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CHAPTER 1 


INTRODUCTION TO AIRCRAFT PARAMETER IDENTIFICATION METHODS 
K, H, Hohenenser and D, Banerjee 


ABSTRACT 

Some of the more important methods are discussed 
that have been used or proposed for aircraft parameter 
identification. The methods are classified into two 
groups: Equation error or regression estimates and 

Bayesian estimates and their derivatives that are based 
on probabilistic concepts. In both of these two groups 
the cost function can be optimized either globally over 
the entire time span of the transient, or sequentially, 
leading to the formulation of optimum filters, Identi- 
fiability problems and the validation of the estimates 
are briefly outlined, and applications to lifting 
rotors are discussed. 
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Preface to Second Yearly Report under Contract HAS2-7613 

Work under Contract NAS2-7613 started on July Ij 
1973 as a continuation of research conducted under 
Contract NAS2-4151 since February 1, 1967, The 
research goals stated in Contract NAS2-7613 are 

(a) Assess analytically the effects of fuselage 
motions on stability and random response. The 
problem is to develop an adequate but not overly 
complex flight dynamics analytics,! model and to 
study the effects of structural and electronic 
feedback, particularly for hingeless rotors. 

(b) Study by computer and hardware experiments the 
feasibility of adequate perturbation models from 
non-linear trim conditions. The problem is to 
extract an adequate linear perturbation model for 
the purpose of stability and random motion studies. 
The extraction is to be performed on the basis of 
transient responses obtained either by computed 
time histories or by model tests, 

(c) Extend the experimental methods to assess rotor 
wake-blade interactions by using a 4-bladed rotor 
model with the capability of progressing and 
regressing blade pitch excitation (cyclic pitch 
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stirring), by using a 4-bladed rotor model with 
hub tilt stirring, and by testing rotor models 
in sinusoidal up or side flow. 

Work on research goal (a) has been reported in 
Part I of the First Yearly Report under subject contract 
titled "Methods Studies' Toward Simplified Rotor-Body 
Dynamics", The results were published in the paper; 
Hohenemser, K, H, and Yin, S, K, , "On the Use of First 
Order Rotor Dynamics in Multiblade Coordinates" presented 
at the 30th Annual National Forum of the American 
Helicopter Society, May 1974, Preprint 831, 

Initial work on research goal (b) has been reported 
in Part II of the First Yearly Report under subject 
contract titled "Computer Experiments in Preparation 
of System '-Identification from Transient Rotor Model 
Tests". 

Initial work on research goal (c) has been 
reported in Part III of the First Yearly Report under 
subject contract titled "Experiments with a Four— Bladed 
Cyclic Pitch Stirring Model Rotor." 

The second Yearly Report under Contract NAS2-7613 
is subdivided into two parts, whereby Parts I and II 
are related to the research goals (b) and (c) 
respectively. The authors and titles of the two 


parts are: 
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Part I Hohenemser, K. H., Banerjee, D. and Yin, S, K. , 
“Methods Studies on System Identification from Transient 
Rotor Tests,” 

Part II Hohenemser, K, H, and Crews, S. T., "Additional 

Experiments with a Four-Bladed Cyclic Pitch Stirring 

\ 

Model Rotor." 

Part I begins with an introduction to aircraft/iden- 
tification methods, it contains the results of computer 
experiments with several selected methods of system 
identification applied to lifting rotors, and contains 
the development of a new method for optimal data 
utilization. 

Part II presents extended frequency response 
measurements with the four-bladed model rotor including 
dynamic wake measurements at zero advance ratio. It 
describes the modifications of the model and its 
instrumentation for transient pitch stirring tests, 
and it discusses the development of software for the 
data processing. 
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Scope of Part I of Second Yearly Report 


The computer experiments presented in the First 
Yearly Report under Contract NAS2-7613 assumed a rotor 
condition with « 4 advance ratio. The input transients 
consisted of rectangular and wave shaped normal flow 
impulses. Noise polluted blade flapping responses 
represented the system output, A linear identification 
scheme was applied that worked well if data preprocessing 
by a digital filter and by a Kalman filter was adopted. 
Blade Lock number and collective pitch angle were 
assumed to be not well known and were identified. 

During the past year this work was extended in several 
directions : 

1, The normal flow transients were replaced by constant 
acceleration pitch stirring transients, since the 
wind tunnel model is particularly suited to accept 
such transients, 

2, For zero advance ratio a two parameter dynamic wake 
model was adopted that included a time constant, A 
single blade representation is now no more applicable. 
A multiblade identification analysis was developed, 

3, The previously used identification method was modified 
and extended to an iterative equation error esti- 
mation with updated Kalman filter. 


REPEODUCIBILITY OF 1.' 
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4, The maximum likelihood method was applied to both 
the single blade and the multiblade identif ication» 

5, The maximum likelihood scheme, assuming the absence 
of system noise, was used to obtain optimum data 
utilization, that is to avoid inaccuracies of the 
parameter estimates from insufficient data pro- 
cessing, and to avoid excessive data processing for 
a given desired accuracy of the parameter estimates. 

Since parameter identification practice is a still 
developing and somewhat controversial discipline, it 
was believed of interest to provide a brief introduction 
to aircraft parameter identification methods in Chapter 
1 of this report. Chapter 2 describes the extensions 
of the work with respect to items 1 to 4, while Chapter 
3 describes the work performed with respect to item 5, 
The 3 chapters are written as independent self suf- 
ficient reports with separate abstracts and lists of 
references for each. 
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The application of methods for system parameter 
identification from transients is a rapidly growing 
field of study. There is a confusing multiplicity of 
methods and of names for these methods; and contradic- 
tory claims are made as to the efficiency and use- 
fulness of the various methods. Widely varying 
parameter identification methods are being used at 
the various aircraft flight test centers. There also 
is a substantial difference between what is theoretically 
possible and what is practically desirable to perform 
a given i ob with adequate accuracy and with a reasonable 
computer effort. Since a comprehensive review of 
aircraft parameter identification methods is not 
available in the lireraturej it was believed useful 
to provide a brief introduction to this field in Chapter 
I before presenting the applications to lifting rotors in 
Chapters II and III, The review of identification methods 
is by no means complete. Only the most important methods are 
discussed. Only rough outlines are given for the 
various methods. Details of the derivations and of the 
application algorithms are found in the cited literature, 

ELEMENTS OF SYSTEM IDENTIFICATION FROM TRANSIENTS 

System identification is the process of extracting 
numerical values for system parameters and other sub- 
sidiary parameters (process and measurement noise 


REPRODUCIBILITY OF T£/ 
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covariances, bias, initial states, etc.) from the time 
history of control or other inputs and of the resulting 
system responses, A schematic for the measurements is 
shown in Fig. 1, The process of system identification 
involves four steps: 

1, Selection of a suitable input that insures parti- 
cipation of all important modes of the system in 
the transient response, 

2, Selection of sufficiently complete and accurate 
instrumentation to measure the key input and output 
variables . 

3, Selection of a mathematical model that adequately 
represents the actual system characteristics, 

4, Selection of an efficient criterion function and 
estimation algorithm for the identification of the 
unknown system parameters. 

The concept of system identification is illustrated in 
Fig. 2, The design input is fed both to the actual 
system and to its mathematical model that contains the 
unknown parameters. The measured response, polluted by 
measurement noise, is compared with the computed 
response from the mathematical model. The difference 
between these two responses, the response error, is used 
in the parameter estimation technique based on the 
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EXTERNAL DISTURBANCE 



Figure 1 Schematic of Measurements for System Identification 



Figure 2 Illustration of System Identification 
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criterion function and optimizing technique. The 
estimation algorithm may also use apriori information, 
e.g,, initial statistics of the parameters. Here we 
will be mainly concerned with the fourth of the 
previously listed steps, that is with the various 
estimation algorithms. 

The mathematical representation of the system will 
be given in the non-linear case by; 

System equation x(t) = f(x,u,t) + r (t)w(t) (1) 

Initial condition x(t = 0) = Xq 


Measurement 

Equat 

ion 

y(t) = h(x,u,t) + v(t) 

(2) 

where x(t) 


n X 

1 

state vector 


u(t) 


P ^ 

1 

input vector 


w(t ) 


q X 

1 

system noise vector, covariance 

Q 

y(t) 


r X 

1 

output or measurement vector 


v(t) 


r X 

1 

measurement noise vector, 
covariance R 


If the system is 

linear. 

equations 1 and 2 reduce to 


x(t ) = 

F(t) 

X Ct ) 

+ 

G(t) u(t) + r(t) w(t) 

(3) 

y(t) = 

H(t) 

X (t) 

•f 

D(t) u(t) + g + v(t) 

(4) 


where g is the vector of bias errors. 
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CLASSIFICATION OF IDENTIFICATION ALGORITHMS 

The various estimation algorithms can be classified 
into two groups presented in Table 1, The first group 
listed in Table 1 above the double line is based on 
statistical regression and does not admit a probabilistic 
interpretation. The algorithms listed in Table 1 below 
the double line are based on probabilistic interpre~ 
tation. In the equation error estimate no measurement 
noise is modeled, the following 4 methods include both 
measurement and system noise, while in the output error 
estimate no system noise is modeled. The various 
algorithms listed in Table 1 will be discussed in the 
following sections. 

EQUATION ERROR ESTIMATES 

Equation error methods assume a performance 
criterion that minimizes the square of the equation 
error (process noise). They are least squares techniques 
and they require the knowledge of all response variables 
(states) and their derivatives. In the so called least 
squares method the unknown parameters are selected 
such that the integral over the square of the state 
equation error is minimized, see for example reference 1. 
With equation 1 we have the error function (the upper 
integration limit T is the time over which the measurements 


are taken) 




Classification 
of Estimate 


Criterion 
for Estimate 


Solution 
for Estimate 


Equation Error Estimate 

(no measurement noise 
is modeled) 


Bayesian Estimate 


0 Minimizes equation error 
squares and/or integrated 
equation error squares 


0 = E(0/Z) 


Closed form global or 
sequential optimization 
by optimum filter (Method 
function 6(t-ti) or e"S't), 


Quasi-Bays ian Estimate 
with augmented state 


Quasi-Bayesian Estimate 


0»x such that 


MAX f(0jx/z; = 
0 , X 


f(0, x/Z) 


0 such that 


Sequential optimization by 
optimum non-linear filter 
(extended Kalman filter 
with or without local 
iteration and/or smoothing) 


Maximum Likelihood Estimate 


MAX f(0/Z) = f(G/Z) 
0 


0 such that 


Iterative global optimi- 
zation with state Kalman 
filter equations as 
constraints e 


MAX f(Z/0) = f(Z/0) 
0 


Output Error Estimate 

(no system noise 
is modeled) 


0 minimizes output error 
squares . , 


Iterative global optimi- 
zation with system 
equations as constraints. 


Table 1 Classification of Estimation Algorithms 


REPRODUCIBILITY OF THE 
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J 


T 



f (x,u,0,t ) 3*^ WEx“f (x,u, 0,t ) 3dt 


(5) 


where 0 is the m x 1 vector of unknown parameters. W is 
a positive definite weighting matrix. An appropriate 

I 

choice for W would be 9 **^ where Q is the covariance of 
the process noise. For the usual digital data pro- 
cessing, the variables x, x, u 'are sampled, and only 
available at discrete time points t^. Mathematically 
the sampling process can be expressed by multiplying 
the system equation with the delta function 6(t - tj^). 

The integral of equation 5 then becomes a sum. One 
can use instead of the delta function also a different 
"method function" , for example e , that would allow 
taking the Laplace transforms. 

If the system is linear in the unknown parameters 
0, the system equation can be written in the form 

x(t) = F(x,u*t)0 t r(t) w(t) (6) 

and the performance criterion (5) becomes 
T 

J = f [x - F(x,u,t )©]'*' W[x - F(x,u,t)0]dt (7) 

o 

Since the function inside the integral has continuous 
derivatives with respect to 0 we set 


3J/90 = 0 


( 8 ) 
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thus resulting in the closed form solution 


0 = 


T T 

^ J' F(x,u,t )dtj”^ J" F"^(x,u,t)W x(t)dt 


(9) 


The first factor is the covariance matrix of the 
estimate. If the system is non-linear in the unknown 
parameters, the solution equation 9 can be replaced by 
an interative solution where F(x,u,t) is substituted 

A 

by 3f (x,u,0k,t )/80 and 9 on the left hand side is 

A 

replaced by shown (for example 

reference 1) that the parameters in the row of 

f(x,u,0,t) are independent of all the elements of x(t) 

o 

except Xj^(t). This Independence is one of the 
drawbacks of the least squares method, in that only 
one of the measured state derivatives is used in 
determining a given rovr of the f(x,u,0,t) matrix. If 
one of the signals has not been measured, the least 
squares method does not provide an estimate of the 
parameters related to that signal. This independence 
also illustrates the fact that the estimate of one row 
of the f(x,u,0,t) matrix is obtained independently of 
the other rows, and no "trade-off” can be made between 
elements in different rows to improve the estimate. 

For some applications it is practical to include 
the state vectors in the error minimization. In the 
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modified least squares method a combination of the 
standard least squares with the integrated least 
squares is used. The parameters obtained by this 
method not only trace the derivative of the state 
but also the state itself over the selected time 
interval. The performance criterion now includes 
in addition to the equation error also the integrated 
equation error: 


J 



t 

F(x,u,t)0 + J' x(r)dT - 



F(x,u, T )0 dx [ i ^dt 


( 10 ) 


where W is a positive definite weighting matrix and where 

T 


I 1 A 1 1 A A W A 


Minimizing the expression, equation 10 results in the 
estimate 


( 11 ) 


i 

/ 


{F(x,u,t) + 

o 'o 


X 

J F(x,u 


,x) dx} W{F(x,u,x) + 


X 

^ F(x,u,t) dx} 


dt 


-1 


( 12 ) 


- T t t 

J {F(x,u,t) + f F (x,u,x )dx}'^W{x(t ) + f x(x) dx}dt 

*0*0 *o . 


1?T5PR0DUCBlLrrY OF the 
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This method has the same row independence of f(x,u,0,t) 
as the standard least squares method. 

Since these methods do not allow for measurement 
errors, they result in biased estimates when this type 
of error does exist. When measurement errors are small, 
as is increasingly the case in modern instrumentation, 
the equation error method is preferrable over other 
methods because of its simplicity. It is widely used 
also when measurement errors are substantial and then 
serves as start-up technique for the output error and 
other iterative methods. 

In many applications, measurements of some of the 
responses or their derivatives are not available. If 
the response, but not the rate of response is 
measured, it is tempting to differentiate the measured 
response. However, the differentiation of measured 
data introduces additional uncertainty so that this 
technique is usually inaccurate. If e is used as 
a methods function, Laplace transforms can be used. 

The estimation then reduces to an algebraic manipulation 
of the data that avoids their differentiation. The 
Laplace transform technique as a substitute of dif- 
ferentiating measurement data is discussed in reference 2, 
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BAYESIAN AND QUAS 1-BAYESI AN PARAMETER ESTIMATES 

In the preceding methods we specified a cost 
criterion J that represented the "loss'* resulting 
from an incorrect estimation of the unknown system 
parameters. The parameters were then selected in such 
a way as to minimize the loss. If a priori probabilities 
exist not only for the measurement errors but also for 
the unknown parameter vector 0 then one can define an 
expected loss and select the parameter vector in such a 
way as to minimize this expected loss. Such an estimate 
is called a Bayesian estimate, see for example reference 3, 
The form of a Bayesian estimate depends on the form 
of both the loss function and of the a priori distri- 
bution of the measurement and parameter noise. For the 
particular case of positive semi-definite quadratic loss 
functions the Bayesian estimate is the mean of 0 
conditioned on the observations. This is true regardless 
of the distribution of measurement and parameter noise. 

It has also been shown that for the case of unimodal 
symmetric a posteriori distribution of the parameters 
given the observations, the Bayesian estimate is the 
conditional mean for a considerably wider class of loss 
functions. For these reasons the Bayesian estimate can 
be defined as the conditional mean of the parameter 


distribution 
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In order to compute the conditional mean it is 
first necessary to determine the conditional probability 
density for 9, This density can be written from Bayes 
rule as (Z is the set of all observations) 

f(0/Z) = f(Z/0) f(0)/f(Z) 

The denominator is a normalizing factor determined from 

f(Z) = / f(Z/0) f(0) d0 

all 0 

The optimal Bayesian estimate is now given by 

0 = (l/f(Z) ) Jb f(Z/0) f(0) d0 
all 0 

In general the evaluation of equations 14 and 15 would 
require the solution of the system equations for all 
possible values of the parameter vector 0. This is a 
large effort, especially if the dimension of 0 is large. 
If f(0/Z) is unimodal and symmetric about its mean 
value» the conditional mean corressponds to the mode. 
Since f(Z) is merely a scale factor the finding of the 
mode requires neither the evaluation of the integral in 
equation 14 nor that in equation 15. The mode 0 of 0 
has the property 

max fiO/Z) ~ max f(Z/0) f(0) = f(Z/0) £(0) 

0 


(13) 


(14) 


(15) 


( 16 ) 



13 


Even if the a priori density f(0) is symmetric it does 
not follow that the conditional density f(0/Z) is also 
symmetric since in general the observations depend 
non-linearly on the parameters. Estimation according 
to equation 16 is, therefore, called "quasi-Bayesian” 
estimation. Another designation used for example in 
reference 4 is maximum a posteriori probability (MAP) 
parameter estimate. Since the logarithm is a monotonic 
function of its argument we can replace equation 16 
by maximizing the expression 

J = log f(Z/0) + log f(0) (17) 

If no a priori information about the parameters 0 
is available, that is, if the a priori density is uniform, 
f(0) = constant, the quas i-Bayes ian estimate reduces to 
the ’’maximum likelihood” estimate which involves 
finding the maximum of f(Z/0), 

ESTIMATES ASSUMING GAUSSIAN DISTRIBUTIONS 

The evaluation of equation 17 becomes particularly 
convenient if we assume Gaussian densities for the 
parameters, for the observations and for the system 
states. In linear systems and linear measurement 
equations (equations 3 and 4) one needs only to assume 
that the system noise w(t) and the measurement noise v(t) 
is Gaussian, It then follows that states x(t) and 
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observations y(t) are also Gaussian, For non-linear 
systems with Gaussian noise f(Z/0) tends to a Gaussian 
density as the sampling rate is increased (see for 
example reference 5 p, 29), The assumption of Gaussian 
densities for all variables is, therefore, a reasonable 
one. Since 0 is a m x 1 vector we now have the a priori 
density 


f(0) 


Pq exp - (l/2)(0-0)’^ (0-0) (18) 


Pg is the parameter covariance matrix, 0 is the 
parameter mean. Except for a constant additive term, 
log f(0) is now given by 


log f(0) = (-1/2) (0-0)^ P“^ (0-0) 


(19) 


In order to obtain an expression for log f(Z/0) in 
equation (17), we assume that Z consists of N consecutive 
observations y(l) .. y(N), 

Z = Yj^ = {y(i), . . . y(N)> (20) 

With successive application of Bayes rule we obtain 
f(Yj^/0) = f(y(l), . . ,, y(N)/0) = F(y (N)/Y^j_^,0) f(Yj^__^/0) 


N 

= IT f (y ( j )/Y . 
j=l ’ 


- 1 ^ 


0 ) 


( 21 ) 


^epeoducibility of the 
SSal page is pooe 
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Taking the logarithm we have 


N 


log f(Yj^/0) = Z log f(y(j)/Y. 0) 

j=l 


( 22 ) 


^y(j)/Yj_i,e) 


is the observation estimate at time j 
given all preceding observations and given the parameters. 
We denote the observation by y( j ) and its expected value 

A 

and covariance respectively by y(j/j-D and B(j/j-l), 

We further denote the ''innovation” by 


yC j ) - y( j / j“i) = v( j ) 


(23) 


Since y(j ) is a r x 1 observation vector, its Gaussian 
density is 

f(y(j)) = |b( j/j-i) exp|-(l/2)v^(j ) 


B“^(j/j-l)v(j)| 


(24) 


Taking the logarithm of equation 24 , summing according 

to equation 22 , inserting in equation 17 and inverting 

the sign we have now to minimize the expression (see 

also equation 19) 

N 
E 

j=l 


S {v"^(j )B“^(j/j-l)v(j ) + log [ B (^‘ / j-1) I > 


+ (0-0)'^ P“^ (0-0) 


(25) 
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If no a priori information is available before taking 
observations, the last term in the expression 25 is 
constant and we then have the criterion for the 
maximum likelihood estimation, Bayesian or quasi- 
Bayesian estimation is rarely used since a priori 
densities for the parameters are in most applications 
not available. 

MAXIMUM LIKELIHOOD PARAMETER ESTIMATION 

According to expression 25 maximum likelihood 
estimation is equivalent to minimizing the so-called 
likelihood function 

J(0) = E [v (j)B"-^(j/j-l)v(j) + log |B(j/j-l)|] 

i=l 

In the presence of system noise the minimization of the 
expression 26 is very difficult. When going from time 
j-1 to time j one first has to solve the prediction 
equations for the estimate of the state and for its 
covariance. Assuming the linear system equation 3 with 
zero mean Gaussian system noise w(t) the prediction is 
given by 

x(j/i-l) = F x(j/j-l) + G u(t) , (j-1) < t < j 
P(j/j-l) = F P(i/j-l) + P(j/j-l)F^ + r Q 


(26) 


(27) 


( 28 ) 
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where Q is the system noise covariance and P the state 
covariance. These equations use the estimated state 

A 

and its covariance at time j-1: x(j-l/j-l) and 

to predict the state and its covariance at 
time j: x(j/-i-l) and P(j/j~l). This is the prediction 

before we know the result of the observations at time 
j. After the observations y(j ) have been made the 
optimum estimate is given by the Kalman filter 
equations for the state and for its covariance: 

x(j/-]) = x(j/j-l) + K(j) (y(j) - H x(j/j-D) 

P(j/i) = (I - K(j) H) P(j/j-l) 
with the filter gain 

K(j) = P(j/j-l) (H P(j/j~l) h'^ + R)“^ 

where R is the measurement noise covariance. The 
covariance of the observations B(j/j-l) that occurs in 
the cost functional 26 is given in terms of the state 
covariance before observations by 

B(j/j-l) = H P(j/i-l) + R 


(29 ) 

(30) 

(31) 

(32) 


Thus the terms in the expression 26 that is to be 
minimized require the solution of the prediction equations 
27 and 28 for each time interval and of the up-date 
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equations 29 and 30 at each sampling time together with 
the solution of the measurement equation 4, Reference 
5 gives an algorithm for the solution of the problem. 
However, due to its complexity this algorithm has not 
as yet been applied to a practical problem of aircraft 
parameter estimation, see for example reference 6, 

The problem of minimizing the expression 26 is 
greatly simplified if the observation covariance 

can be assumed constant. This is for example 
true for zero system noise, when according to 
equation 32 B(j/j-l) = R. The problem then reduces 
to minimizing the cost function 

^ T -1 

J(0 ) E V ( j ) R v(i ) 
j=l 

where v(j) is given by the innovation term 23, Since 
equation 33 represents (according to equations 2 and 4) 
the sum of the measurement error squares, the estimation 
with equation 33 is also called output error method of 
estimation. There are several algorithms available to 
perform the optimization of J(0) from equation 33, The 
most widely used is the modified Newton-Raphson or 
quasiliniearization method. It has the advantage that 


(33) 


the sensitivity or information matrix is obtained as a 
byproduct. The inverted information matrix gives the 
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Cramer-Rao lower bound for the parameter covariance. 

This lower bound is found in many applications to be a 
more meaningful measure of the accuracy of the parameter 
estimate than the parameter covariance obtained from 
the equation error method (first factor in equation 9). 

OUTPUT ERROR METHOD USING QUASILINEARIZATION 

We use an iterative method beginning with an 

A 

initial parameter estimate 0 = 0 q . The problem is to 
find a zero of the gradient of the cost function 33, 

3J/30 = 0, Consider a two-term Taylor series 
expansion of 3J/30 about the kth value of 0 

(3J/30)3^+1 ~ ^ O^J/90^)k A 0]^+jL 

where . 

A 

^ ®k+l = ®k+l ~ ®k 

2 2 

(a J/aO )j^ is the second gradient of the cost function 
with respect to 0 at the kth iteration. If equation 34 
is a sufficiently close approximation, the change in 0 
for the ^k + ijth iteration to make (aJ/30)]^+2 approximately 
zero is 

A 0}^+i = - C(3^J/30^)j^]“^ (aj/a9)k (35) 

Using now for v(j) the two term expansion 

v(j)k = ^(i)k-i + ^ ®k 


(36) 
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one obtains for the first and second gradients of the 
cost function 

(9J/30)j, = 2 E ^ R Cv(j)]i, 

i = i 


M 






g 

We thus need solutions for v(j)j^ and — v(j)j- . For 
this purpose we first solve the system’ and measurement 
equations 

X = f(x, u,t) 

A ^ 

y = h(x,u,t ) 


for each iteration whereby the initial conditions are 
either obtained from the measurements or are included 
in the unknown parameters 0. The innovation is now 
obtained from equation 23. Next we solve the 
"sensitivity equations" for each iteration 


3x/30£ 


3f/30. + 3f/8x 


A 

x = x 


3x/30£ 


3y/30£ 


3h/3x ^ 3x/30^. 

x=x 1 
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( 37 ) 


(38) 


(39) 


(40) 
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The initial conditions of 3x/30^ are zero except when 
x(0) is identified as part of the parameters 0. In 
this case the initial partials have the value one. With 
equation 23 we can now compute the first and second 
gradient of the cost function* equations 37 and 38, 
and then obtain the change in parameters for the next 
iteration from equation 35 . This involves the 
inversion of the sensitivity matrix M (equation 38), 
whereby is the Cram^r-Rao lower bound for the 

covariance of the parameters. 

The method is easily extended to the case with 
a priori information on the parameters, equation 25, 

The sensitivity matrix 38 is then augmented by the 
term 2Pg^, and the gradient 37 is augmented by the 
term 2Pg^ (Oq - see reference 2, 

PARAMETER ESTIMATION BY FILTERING 

The parameter estimation methods discussed so far 
can be denoted as “global" methods. The performance 
criterion includes the test data for the entire duration 
of the transient. Filtering is an important tool in 
state and parameter estimation. It can be used either 
in conjunction with global estimates, or it can be 
used as a direct approach to state and parameter 
estimation. An example of the first type of filter 
applications is the prefiltering of test data before 
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using them in a least squares regression estimate, 
see for example reference 4, The Graham digital filter 
can remove high frequency noise. A Kalman filter 
can be used to estimate state variables and their 
rates not directly measured. It also removes the 
noise in the measurements. The role of the Kalman 
filter in maximum likelihood estimation has been 
shown in equations 29 to 31, where it is used to 
establish the innovation sequence. 

In addition to applications in global estimation 
methods, filters can also be used as substitutes for 
global methods. The advantage of such direct filter 
methods is a reduction in computer effort particularly 
in cases with a large number of parameters. The dis- 
advantage is that unlike the inverted information 
matrix of the maximum likelihood method that provides a 
lower bound on the parameter covariances, no physically 
meaningful parameter covariances are obtained with the 
direct filter methods. The covariance propagation 
equations require initial values that are usually 
impossible to obtain in any rational way. Though 
improvements of the filter solution (forward time 
integration) can be achieved by smoothing (backward 
time integration), the final parameter covariances 
remain arbitrary, since they evolve from arbitrary 
initial covariance estimates. 
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Assuming that all state variables and their rates 
have been either measured or are otherwise known from 
manipulating the measurement data, the unknown 
parameters, if they occur in linear form in the 
state equation, can be found by application of a 
linear filter, see for example reference 7. The 
classical regression method is a special case of this 
direct filtering method, namely for infinite initial 
parameter covariances* In classical regression one 
obtains a single value of the error covariance matrix. 
The direct filter application allows the use of a 
finite initial error covariance matrix and it gives 
the evolution of this matrix as a function of time. One 
thus obtains an indication when to stop processing the 
test data after their information contents has been 
exhausted. As mentioned before, the absolute values 
of the error covariances are meaningless, since one 
usually does not have a rational way of establishing 
initial values for the parameter covariances, 

A method that appears to be economic of computer 
time for large numbers of unknown parameters was used 
in reference 4 for application to helicopters. The 
method consists of a simultaneous identification of 
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states and parameters with the help of a non-linear 
filter. In other words, the unknown parameters are 
treated as additional state variables. Since there 
occur products of state variables and parameters, the 
filter is a non-linear one. The so called extended 
Kalman filter appears to be particularly useful for 
this purpose. Either non-linear filtering alone or 
linear filtering in combination with smoothing is 
performed. The absolute values of the parameter 
covariances are again of no physical significance since 
they depend on the arbitrary initial values. In the 
following, a brief discussion of the direct use of 
filters in parameter estimation is given, 

LINEAR FILTER METHOD OF PARAMETER ESTIMATION 

The equation error estimate based on the perfor- 
mance criterion of equation 5 was given in global form 
by equation 9, which is valid provided that the system 
is linear in the unknown parameters 0, The global 
estimate can be replaced by a filter solution. 

Consider a system with the process equation 

0 = 0 

and with the measurement equation 

x(t) = f(x,u,e,t) + v(t) 


(41) 


(42) 


REPRODUCIBILITY OF THE 
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The extended Kalman filter that minimizes sequentially 
the performance criterion 5 is given by the filter 
equation 


X 

0 = 


PgOf/90) 


W (x - (9f/3O)0) 


(43) 


and by the covariance equation 

Pg = -PgOf/SO)’^ W Of/30) Pg 


(44) 


x and X are here quantities that are measured or that 
have been reconstituted from measurements by prefiltering. 

The weighting matrix W can be interpreted as the inverse 
of the error covariance for v(t). The sequential 

estimate from the filter solution- gives the optimum 
based on all preceding data. One can either select 
an initial estimate 0(0) to integrate equation 43, 
or one can also assume a zero value for the initial 
estimate. The initial error covariance Pg(0), as 
mentioned before, is usually not known. One should, 
therefore, use a rather large value for the initial 
covariance. Though the resulting values of Pg(t) obtained as 
function of time by solving equation 44 have no physical 
meaning as far as absolute values are concerned, Pg(t) 
will asymptotically approach zero, and it can be used as 
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a criterion when to stop processing the data. As can be 
seen from equation 43, no further changes in the 
parameters 0 will occur, once Pg has approached zero. 

The covariance equation 44 has the closed form 
solution (see reference 7) 


t 

‘Pg(t) = [Pg^ (0) + 

o 

and the filter equation 43 has the closed form solution 


/ 


(3f/30)^ W (3f/30) dT] 


-1 


0(t ) 


Pq<« 


/ 


[p“^ (0)0(0) + / (3f/30)" W X dT] 


Whether it is computationally simpler to numerically 
integrate equations 43 and 44 or to evaluate equations 
45 and 46 depends on how many points in time are to be 
covered. For infinite initial parameter covariance, 
Pg^(O) = 0, equation 46 becomes identical to equation 9, 

BAYESIAN ESTIMATION AS A FILTERING PROBLEM 

If we extend the quasi-Bayesian or maximum a 
posteriori probability (MAP) criterion 16 to include 
both the parameters 0 and the states x(t), we have 

max f(x(t),0/Z) - max f(Z/x(t),0) f(x(t),0) 

x(t),0 x(t),0 


(45) 


(46) 


( 47 ) 
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Assuming now the non-linear system and measurement 
equations 1 and 2, and assuming further that states 
and parameters have Gaussian distributions of the 
form of equation 18, the criterion 47 becomes (see 
for example references 4 and 8) one of minimizing 
the quadratic function 


J 


( 1 / 2 ) 



+ 



ly (t )-h(x,u ,t ) 


^ + I |w(t) 

R-l 




subject to the constraint equation 1, If the system and 
measurement equations 1 and 2 are linearized about the 

^ A 

current estimates x and 0 the recursive solution of the 
minimization problem 48 results in the extended Kalman 
filter equations given for the continuous case by (see 
reference 4) 

X = Of/3x)x + Of/3u)u + P(3h/3x)^R“^(y-(3h/3x)x) 

0 = P (3h/3x)’''’R"^ (y - (3h/3x)x) 

p = (3f/3x)P + POf/Sx)”^ - P(3h/3x)'^R“^(3h/3x)P + Q 
+ (3f/30)P0j^ + C(9f/30)pQx^'^ 

P0X = P0j^(3f/3x)'^ + P»(3f/30)^ - p0xOh/3x)^R”l(3h/3x)P 

P’ = - P0xOh/3x)'^ R'^(3h/3x)P^^ 


(48) 


( 49 ) 


with 
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P is as before the state covariance, P 0 the parameter 
covariance, and the combined covariance where 

the parameters 0 are included in an augmented state, 

P’ is an abbreviation given in the last of equations 
49, Even if the original system is linear, the 
augmented system is non-linear and hence the filtering 
problem must be solved by a non-linear filtering 
technique. In reference 4 the raw data are preprocessed 
by a digital filter and by a Kalman filter that does 
not use the unknown parameters but merely makes use 
of the transformation equations from a space-fixed to 
a body-fixed reference system (Euler equations), 

bebacqz in reference 9 applies basically the same 
method except for a discrete instead of the continuous 
filter formulation. He further uses a one stage 
filtering-smoothing algorithm which has the advantages 
of reducing the bias due to non-linearities and of making 
the algorithm less sensitive to initial conditions, 

Mehra in reference 4 is critical of using an extended 
Kalman filter for the augmented state including the 
unknown parameters. His arguments are that the 
uncertainties in the states are usually much smaller 
than the uncertainties in the parameters. Therefore 
the assumption of local linearization about the latest 


tBt 
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estimate which are acceptable for state estimation with 
an extended Kalman filter are generally less valid for 
parameter estimation<> Moreover, the filter for the 
augmented state assumes knowledge of the a priori 
parameter covariances which are unknown. As mentioned 
before, the arbitrary a priori parameter covariance 
used as initial conditions for a filter that includes 
parameters as state variables gives unreliable confidence 
limits on the parameter estimates. An added difficulty 
of applying a filter to the augmented state is that poor 
a priori estimates of the parameters make the convergence 
rate slow or may even cause divergence of the filter 
solution. Though improvements can be applied to the 
extended Kalman filter like local smoothing and local 
iteration and smoothing, the basic shortcomings of this 
method appear to have been correctly described in reference 4, 
Unfortunately, the application of the complete algorithm 
of maximum likelihood identification given in reference 
5 is for a large system much more demanding of computer 
size and time than the filter solution with the 
augmented state. While aircraft parameter identification 
with the complete maximum likelihood algorithm of 
reference 5 has not as yet been accomplished, the method 
of filtering the augmented state has been applied to 
several aircraft parameter identification cases, for 
example in references 4 and 9, 
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IDENTIFIABILITY PROBLEMS 

Identif lability problems can occur no matter what 
identification algorithm is used. They are related to 
the initial 3 steps involved in system identification 
as listed at the beginning of this chapter: the 

selection of a suitable input, the selection of the 
instrumentation, and the selection of the mathematical 
model. A few comments are added here to point out some 
difficulties that have been encountered due to these 
three initial steps. 

If the input does not adequately excite some of 
the system modes, the associated parameters cannot be 
adequately identified. Sometimes it is practical to 
combine the responses to various types of inputs into 
a single identification run, see reference 4, While 
each of the single inputs excites only a limited number 
of modes, the combination of inputs provides an adequate 
excitation of all modes required for the estimation of 
the parameters. Efforts have also been made to design 
inputs on the basis of certain optimization criteria. 

More details on this problem are given in Chapter 3 of 
this report. 

If there are large unaccounted for instrumentation 
errors non-physical parameter values may result. In 
reference 10 instrumentation lags and control measurement 



31 


errors were found to be most significant. Static 
measurement errors and instrumentation lags can be a 
much greater source of narameter inaccuracies than 
white noise. A detailed analysis of the relationship 
between static and dynamic measurement errors in states 
and control inputs and the accuracy of the parameter 
estimates is required. 

If the selected mathematical model for the system 
is inadequate, the parameters are forced to account for 
some unmodeled effects. The estimated parameters may, 
therefore, be quite different from those determined by 
aerodynamic theory or wind tunnel tests would indicate, 

A good example is given in reference 11 where a six 
degree of freedom mathematical model for a helicopter 
gave unrealistic derivatives, since it had to account 
for effects of some neglected modes. A unique six 
degree of freedom linear model for the helicopter flight 
dynamics does actually not exist. When a nine degree of 
freedom mathematical model is used, these difficulties 
disappear. Modeling errors are also a major cause for 
the lack of convergence of iteration procedures or of 
parameter identification by filtering methods. The best 
remedy against difficulties from modeling errors is the 
adoption of a more suitable mathematical model. Some 
other measures to improve the convergence of iteration 
procedures or of filtering methods will be briefly discussed. 
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In the cases where a priori values of parameters, for 
example from theory or from wind tunnel tests, are 
available, one can use an a priori weighting matrix 
that expresses the confidence in these values and 
prevents the algorithm from deviating too much from 
the a priori values. Sometimes there exist some 
relationships between the parameters.^ These should 
then be used as constraints in the optimization problem 
to avoid non-physical parameter estimates. If para- 
meter dependencies exist, difficulties are encountered 

I 

in inverting the information matrix. An exact 
dependency between parameters should result in a zero 
eigenvalue of the information matrix, A rank deficient 
solution makes use of the fact that in case of near 
parameter dependencies there is a large spread between 
a set of small eigenvalues and another set of much larger 
eigenvalues of the information matrix. 

In filter solutions, divergence because of modeling errors 
can occur when the covariance matrix becomes prematurely 
too small, thus preventing further test data to be of 
influence. There are several ways to prevent premature 
small covariances. One can provide fictitious noise 
input to the system or one can directly increase the 
parameter covariance in each time step according to 
some rule. One can also overweigh the most recent data 
thus causing the filter to reduce its memory of the data 
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of the more distant past. This indirectly increases the 
parameter covariance matrix. Since too short data length 
and too large errors in the initial parameter estimates 
may also result in non-physical parameter values or in 
divergence of the identification algorithm, longer 
transients and better a priori parameter estimates can 
lead to the avoidance of these difficulties, 

VALIDATION OF ESTIMATES 

Once a set of parameter estimates has been obtained 
the question arises; what confidence can be associated 
with this set? As mentioned before, the parameter 
covariance matrix obtained by filtering the augmented 
state is not a good measure of this confidence. The 
inverted information matrix obtained with the maximum 
likelihood method represents the Cramer-Rao lower bound 
for the parameter covariances and is a better measure 
of this confidence. Using the parameter estimates to 
predict the transients from which the estimates have 
been obtained, and computing the rms error with respect 
to the measured transients, gives another confidence 
measure. However, if the system is inadequately modeled, 
one may obtain a small rms error despite the fact that 
the parameter values are wrong in comparison to theoretical 
or wind tunnel results, see reference 11, A better way 
of validitation is to compare the prediction with the 
results of test data not used in the identification process. 


In fact, it is good practice not to use all of the 
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available test data for the parameter identification but 
to reserve some of the runs for such a comparison. 
Sometimes it is desirable to perform the parameter 
identification not just with one mathematical model but 
with a variety of models. In the case described in 
reference 11 a mathematical model with more parameters 
gave a much better identification result than a model 
with fewer parameters, better in the sense of an 
improved correlation with theoretically and wind tunnel 
generated parameters. However, there are also cases 
where mathematical models with a larger number of 
parameters gave worse identification results than a 
model with fewer parameters, see reference 12, 

Adequate parameter estimation from transients requires 
careful attention to the many contributing factors in 
the input, instrumentation, mathematical modeling, and 
the estimation algorithm; and the validation of this 
process can only be considered complete after the rms 
errors of the prediction with the estimated parameters 
as compared to test data have been found acceptably 
small for all types of possible transient excitations 
of the system. 


^EPRODUCBILrrY OP TBE 
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APPLICATIONS TO LIFTING ROTORS 

I 

Contrary to the fixed wings of airplanes, lifting 
rotor characteristics are not well approximated by the 
usual set of aerodynamic derivatives. One reason is 
.blade modes that must be considered particularly in 
rapid transients. Another reason is the dynamic rotor 
wake that is produced by the time varying rotor thrust 
and rotor pitching and rolling moments and that has a 
feedback effect on the rotor forces and moments. The 
omission of the blade modes, as shown in reference 11, 
results in non-unique and non-physical, rotorcraft 
derivatives. The identification is better if separate 
rotor degrees of freedom are introduced even in the 
crude form of a first order lag as was done in reference 
13. 

A variety of identification methods have been used 
with respect to lifting rotors. After preprocessing the 
test data with a digital filter followed by a Kalman 
filter that does not contain the aerodynamic derivatives 
(transformation or Euler equations), least squares iden- 
tification is applied to rotorcraft transient flight 
test data in references 4 and 11, Each identification 
run is made with several transients simultaneously. The 
least squares results are then used as start-up values 
for an extended Kalman filter for the augmented state. 

It is not obvious that the extended Kalman filter 
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actually improves on the least squares results, though 
filter convergence is achieved. In reference 13 the 
output error method with quasilinearization is applied 
without preprocessing the flight test data. The flight 
data of both references 4 and 13 were obtained in calm 
air. The equation error method in its filter form was 
applied in reference 7 to simulated noisy blade flapping 
and torsion measurements at high rotor advance ratio. 

The simulated data were preprocessed by a Graham 
digital filter, but not by a Kalman filter. Reference 7 
assumed that all states and their derivatives had been 
measured. In contrast reference 14 assumed that only 
flapping deflections are measured but not flapping rates 
or flappin g ac celerations. For the dynamic wind tunnel 
tests simulated in reference 14 there is no way of 
applying a Kalman filter thax does not contain the 
unknown parameters. However, it was found that for the 
cases studied, a Kalman filter with considerable errors 
in the unknown parameters was useful in obtaining the 
non measured flapping rates and accelerations. The 
parameter identification was then performed by the 
equation error method in its filter form. 

In Chapter 2 of this report the same method (except 
for using global estimates) is used in an iterative form. 
In addition, the output error method with quasilineari- 
zation is applied to the same and to more complex rotor 
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identification problems. As will be shown, the first 
method is more computer cost effective in cases where 
it works, while the second method is more reliable 
and more versatile. 
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Chapter 2 


COMPUTER EXPERIMENTS WITH PITCH STIRRING TRANSIENTS 
K, H, Hohenemser and S. K« Yin 


ABSTRACT 


Selected methods are applied in the form of computer 
experiments to two problems of lifting rotor state 
variable and parameter identification. The first 
problem refers to a rotor condition at .4 advance ratio. 
Cyclic pitch stirring with constant acceleration is 
assumed as known input. Noise polluted blade flapping 
measurements are assumed as the only measured output. 

Blade Lock number and collective pitch are the unknown 
parameters to be identified. The second problem refers 
to a rotor condition at zero advance ratio again with 
constant acceleration pitch sitrring, A two parameter 
dynamic inflow model is stipulated. The only measurements 
are noise polluted blade flapping responses. Blade Lock 
number and the two parameters of the dynamic inflow model 
(Including a time constant) are the unknown parameters 
to be identified. 

Two parameter identification methods are applied and 
their relative merits evaluated: Iterative equation 

error estimation with updated Kalman filter and the 
maximum likelihood method. The latter method, though 
requiring increased computer effort per iteration, was 
found to be more versatile. 
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INTRODUCTION 

The identification method used in reference 1 is 
based on the experience gained from references 2 and 3 
where states and parameters were determined simul- 
taneously with the help of an extended Kalman filter. 
The filter can easily diverge unless good initial 
estimates are available. Particularly in reference 
3 a considerable effort was applied to obtain such 
good initial estimates. The test data were first 
processed with a digital filter that took out high 
frequency noise without distorting the main signals. 

The data were then processed with a Kalman filter based 
on the Euler equations, which do not contain the 
unknown parameters. Thus measurement bias was removed 
and missing channels were reconstituted. Finally a 
least squares algorithm was applied to obtain estimates 
of the unknown parameters. The subsequent application 
of the extended Kalman filter led to modified parameter 
estimates, however it is not clear whether or not 
these modifications represent improvements. In any 
case the modifications were not large, and the initial 
estimates appeared to be satisfactory approximations. 

In trying to apply the experience from reference 3 
to wind tunnel model transients a difficulty arises. 
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in that no equivalent to the Euler equations for the 
aircraft exists. Thus there is no way of using a 
Kalman filter which is free of the unknown parameters. 
Instead, if a Kalman filter is to be applied, 
estimates of the parameters must be inserted. Another 
difficulty for our wind tunnel model tests is that 
only flapping deflection measurements are made, while 
the rates of deflection and the accelerations are not 
measured. Thus the Kalman filter with the estimated 
parameters is called upon to provide both rates and 
accelerations. Finally, the least squares algorithm 
of reference 3 was replaced in reference 1 by a 
sequential linear estimator for the parameters. 

This has the advantage that finite initial parameter 
covariances can be used, and that the time history of 
the parameter covariance provides a measure for the 
time beyond which no more useful information can be 
extracted from the test data. As far as digital 
filtering is concerned it was found in reference 1 that 
a cut-off frequency range from 2.5 to 2.9 removed the 
high frequency noise without unduly distorting the 
main signal. 

The method used in reference 1 worked well in the 
computer experiments using normal flow transients at 
,4 advance ratio for the identification of the Lock 
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number and of the collective pitch angle. However, a 
number of questions had to be answered by the sub- 
sequent studies. It had been decided to begin wind 
tunnel model rotor transient testing with constant 
acceleration pitch stirring inputs, and later follow 
up with normal flow transient testing. Therefore, it 
was desirable to perform computer experiments with 
pitch stirring transients. All of the identification 
work reported here concerns such transients. The 
linear sequential estimator used in reference 1 
requires the simultaneous integration of the filter 
and of the covariance differential equations, A 
simpler ••global” estimate requires only the inversion 
of a system of linear equations for the unknown 
parameters and the evaluation of a number of integrals 
over the time period of the transient. Therefore, a 
number of comparisons were made between these two 
methods. Finally, the Kalman filter for the test data 
requires estimates of the unknown parameters. The 
question arises whether an iterative form of the 
method is practical, where the identified parameters 
are reinserted into the Kalman filter and a second 
identification is performed, etc. This iteration 
method was tried out in several cases. 


WRODUCBILITY OF THE 
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While the substitution of the theoretical blade 
Lock number by an equivalent Lock number can be 
expected to provide a reasonably good approximation of 
dynamic rotor wake effects if the transient is 
relatively slow and does not contain high frequenclesj 
a better dynamic wake representation is given by 
introduction of the L-matrix (reference 4) together 
with a time constant (reference 5), For zero advance 
ratio the L-matrix degenerates into a L-scalar, so that 
a two parameter dynamic wake description is obtained, 

A single blade representation is now no more applicable, 
and a multi-blade analysis is necessary, A number of 
computer experiments were conducted for this case 
whereby the blade Lock number was treated as third 
unknown parameter. 

During the last decade the maximum likelihood 
method of parameter identification has been success- 
fully applied to airplane and helicopter transient testing. 
This method does not require preprocessing of the test 
data and also does not need complete measurements of 
the deflecitons, of their rates and of the accelerations. 
The parameter covariance estimates obtained with this 
method are more meaningful than those obtained with 
the linear sequential estimator used in reference 1, A 
number of cases were treated with the maximum likelihood 


method 
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In its most general form presented in reference 6 
the maximum likelihood method is capable of handling 
identification problems for cases with both measurement 
and system noise. The method then becomes computa- 
tionally quite complex and has not as yet been applied 
to aircraft transient testing. The usual assumption 
is that the system is free of noise. In this case 
the maximum likelihood method is greatly simplified 
(output error method). It has been applied in this 
report in the simplified form. If the random system 
noise is measured, the simplified method is still 
applicable. However in the planned wind tunnel model 
testing the system noise will not be measured, and 
though it is expected to be small, the question is 
whether or not it will unacceptably degrade the 
estimates of the parameters. Therefore a few computer 
experiments were conducted with both simulated mea- 
surement noise and with simulated system noise, whereby 
the maximum likelihood method was applied in its 
usual simplified form, 

EXCITATION OF PITCH STIRRING TRANSIENTS 

For the wind tunnel experiments with pitch stirring 
transients the initial state of the rotor will be given 
by prescribing the advance ratio, the collective pitch 
angle, the rotor angle of attack (set at approximately 
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zero) and the cyclic control setting that will be 
zero longitudinal cyclic and 1.5® lateral cyclic* 

At the time t^ pitch stirring is initiated. Denoting 
by ( 1 ) the angular pitch stirring speed, positive in 
the direction of rotor rotation, and by m the pitch 
stirring angular acceleration assumed to be constant, 
we have 

w = 0)(t - t^) (1) 

For a progressing mode m is negative, for a regressing 
mode u is positive. In a rotating reference system the 
blade pitch angle is given by 


G = 0Q + 1.5 cos Cu{t-tQ) t t] 

' 

0 for t ^ tg 

“ ■ , ( 2 ) 
oj(t - to) for t > to 

In a multiblade representation the blade pitch angle of 
the kth blade is 

0]^ = 0Q - sin + 0JJ cos (3) 


where 


01 = 


1.5 sin uCt-tg) 


1.5 


'II 


1.5 cos U3(t-tg) 


for t 4 t^ 
for t > to 
for t 4 t^ 
for t > tQ 


(4) 


(5) 
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0j represents forward cyclic pitch, represents 

left cyclic pitch. The wind tunnel experiments will be 
conducted with a variety of pitch stirring accelerations. 
The computer experiments were mostly conducted with a 
pitch stirring acceleration of 

w = - ,l/ir (6) 

which is in the progressing sense. 

Since in the non-dimensional time units used here 
the time of one rotor revolution is 2ir , the angular 
pitch stirring velocity one rotor revolution after 
initiation of pitch stirring is ,2, that is one fifth 
of the rotor angular speed. Figure 1 shows the time 
history of blade pitch for about two rotor revolutions 
(tQ = 12, t = 12 to 24) in a rotating frame of reference. 
Figure 2 shows the time history of blade pitch in 
multiblade representation, that is 0^ and 0^.^ vs, 
time t , Figures 1 and 2 refer to the progressing mode, 

A convenient way of conducting the computer 
experiments is to impose at time t = 0 a step input 
of lateral cyclic control. If to is sufficiently large, 
the transient from the step control input will have 
subsided when pitch stirring begins at t = to* It 
was found that for single blade identification a value 
of to = 12 is adequate. In the multiblade identification 
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including the time delayed rotor wake this value was 
found to be insufficient. Instead, tg = 70 was used. 
Measurement noise was simulated by polluting the 
analytical flapping response with zero mean white 
Gaussian computer generated random sequences. System 
noise was simulated by adding to the multiblade inflow 
coordinates in the system equations zero mean white 
Gaussian sequences. Data processing extended over 
nearly two rotor revolutions, from t^ to t^ + 12, In 
Chapter 3 it will be shown that this choice leads to 
approximately optimal data utilization. If a shorter 
time period were used, large errors in the parameter 
estimates would occur. If a longer time period were 
used, the additional data processing would be unnecessary 
in view of the adequate accuracy of the parameter es- 
timate using a time period of two rotor revolutions. 

For the wind tunnel experiments a variety of collective 
pitch settings 0^ will be tested. For the computer 
experiments the collective pitch setting was mostly 
assumed to be Gq = 2° , 

LINEAR SEQUENTIAL AND GLOBAL FSTTMATORS 

In reference 7 the parameter identification is 
performed from a "system equation" 


0 = 0 


(7) 



and a "measurement equation" 


? = h(5,e) + V 

Equation 8 is actually the system equation arranged 
in a form where the left hand side contains all terms 
that are free of the unknown parameters 0, If the 
system equation is linear in the state variables C and 
in the unknown parameters 0, h(C»0) is a linear 
function of the parameters. The noise vector v refers 
only to the terms on the left hand side of equation 8, 
The state variables that are multiplied by the unknown 
parameters in h(5,0) must be noise free. To obtain 
the 0 both Z and g must be known. If only part of the 
variables in Z and 5 have been measured, Kalman 
filtering is required in order to reconstitute the 
missing terms. 

Denoting by Pq the parameter covariance matrix 
and by R the noise covariance matrix, assuming that v 
is zero mean Gaussian white noise, optimal parameter 

A 

estimates 0 can be obtained by minimizing the cost 


function 
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where theapriori estimates G(0), Pg{0) are assumed to 
be given together with the noise covariance matrix R, 

The differential equations associated with this 
optimal problem are (see for example reference 8) 

0 = P^Oh/GG)'^ R“^(C - h(5,0)) (10) 

pQ = - P.Uh/GG)’^ R”^Oh/30)P-, (11) 

Beginning with the initial a priori estimate for the 
parameters 0(0) and their covariance matrix Pg(0), 
these equations can be integrated and result at each 
time t in the optimal parameter estimate given the 
preceding measurements. Since the initial parameter 
covariance is usually not known and the assumed 
values are rather arbitrary, the matrix P^ from 
integrating equation 11 is not a useful measure of the 
actual parameter covariance. However, once P^ has 
approached zero, the effect of any further measurements 

A 

on the estimate 0 also approaches zero as is evident 
from equation 10, Pg , therefore, is valuable in judging 
for what length of time the data should be processed. 
Equations 10 and 11 represent the "linear estimator" 
used in references 1 and 7, 
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Instead of the sequential estimation by integrating 
equations 10 and 11 with some initial estimates 0(0) 
and Pg(0), one can also obtain a "global" estimate 
directly from equation 9, If one assumes that one and 
the same parameter estimate 0 is valid throughout the 
time range from 0 to T, one obtains by setting 
3J/30 = 0 



Oh/3G)'^R"^(ah/9 0)dt]"^ 


i;Pg^O)0(O) 


+ 



Oh/80)'^R"^ 5 dt] 


( 12 ) 


see for example the appendix of reference 7, A con- 
venient assumption is P”^(0) = 0, which means an 
infinite initial parameter covariance matrix. The 
initial estimate 0(0) is then not required and the 
evaluation of equation 12 is reduced to the deter- 
mination of fixed boundary integrals, a matrix 
inversion and a matrix multiplication. The parameter 
covariance matrix at the time T is given by the first 
factor of equation 12; 

T 

Pq(T) = [Pq^‘ 0) + y Oh/30)'^R“^Oh/90)dt]‘‘^ 

o 


( 13 ) 
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which follows from the integration of equation 11, (see 
for example the appendix of reference 7), P^CT) from 
equation 13 can again be used to judge whether or not 
all the significant information contents has been 
extracted from the data, A comparison between the 
estimation with equations 10 and 11 and with the 
"global'* method of equation 12 will be given later 
for a specific example, 

ITERATIVE EQUATION ERROR ESTIMATION WITH UPDATED 
KALMAN FILTER 

When using the parameter estimation methods of the 
preceding section it is necessary to first determine 
from the noisy deflection measurements estimates for the 
deflections, for their rates, and for the accelerations. 
In reference 1 this was done by passing the noisy 
deflection data through a digital filter that takes out 
the noise above a certain frequency without distorting 
the signal in the low frequency range. The filtered 
deflections were then either differentiated twice, or 
a Kalman filter was applied in order to obtain the 
derivatives. Later studies showed large errors in the 
parameters for too low cut-off frequency of the digital 
filter. It was then decided to omit the digital filter 
and instead use the Kalman filter in an iterative way. 
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In typical examples already the second Iteration was 

as accurate as the result with the combined digital and 

Kalman filter given in reference 1, 

Only simulated noisy blade flapping measurements 

were used in the Kalman filter. The filter provided 

the deflection rates and accelerations needed for the 

’’global" parameter estimate, but not the deflections 

themselves. In other words, the parameter estimate 

was performed with the simulated noisy deflection 

measurements and with the rates and accelerations 

from the Kalman filter. In the first iteration a 

Kalman filter with estimated parameter values was 

used (typically 20% error). After updated parameter 

values had been obtained, a second pass with an 

updated Kalman filter was performed, etc. The deflection 

data remained the same for each iteration, but the 

rates of deflection and the accelerations were updated. 

As will be seen in the numerical examples this method 

worked well for the single blade identification. The 

reason for the good parameter identification in spite of 

substantial measurement noise is probably the following. 

} 

In the equation of single blade flapping the flapping 
deflection term includes a large part that is 
independent of the Lock number, and only a relatively 
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small part that is dependent on the Lock number. 
Therefore, the major part of the noisy deflection 
term can be included in the left hand side of equation 
8 in 5 , The method allows noise in this part. The 
smaller part of the flapping deflection term that is 
of aerodynamic origin and includes the blade Lock 
number is part of the right hand side of equation 8 
and should be noise free. Apparently the actual 
noise in this relatively small term does not lead to 
a substantial bias in the parameter estimate, 

MAXIMUM LIKELIHOOD ALGORITHM 

The maximum likelihood method for our particular 
case pertains to a system equation (zero system noise) 

X = f(x,u,0) 

where x is the state vector, u the input vector as- 
sumed to be known. 0 is the vector of unknown para- 
meters that may include initial values of state 
variables. The measurement equation is assumed to be 
linear and of the form 

y = H X + V 

y is the vector of observed quantities, H is a matrix 
relating the state variables to the observations, v 


( 14 ) 


( 15 ) 


is the vector of random measurement errors, assumed to 
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be zero mean white noise with given covariance matrix R 

E[v(t)v'^ (t)] = R 5(t, t) (16) 

R is assumed to be constant with time. Though the 
preceding equations do not show bias terms, bias 
errors could easily be included in the unknown 
parameter vector 0. 

A sample of measurements yj_ y2 * * * 
made during the time of the transient and the parameter 

A 

estimate 0 is selected such that the conditional 
probability of this sample of measurements given 0 
is maximized, 

A 

0 = max p(y^ , . , y^j/®) (17) 

The following steps lead to the maximum of the likelihood 
function p(yj^ * , , y^/Q)» though there is no assurance 
that the maximum is global. The method outlined here is 
called quasilinearization with the modified Newton- 
Raphson method. It assumes Gausian distributions of 
the random variables, 

A 

1 Select an initial parameter estimate 8 = 0Q* 

2 Solve the system equation 14 with this parameter 
estimate 

AAA 

X = f (x , u, 0 ) 


( 18 ) 
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The initial conditions can either be obtained 
from the measurements, or, where this is not 
feasible, they can be included in the unknown 

A 

parameter vector 0. 

3 Calculate for each measurement the "innovation term" 

v.=y.-Hx. (19) 

3 3 3 

4 Solve the "sensitivity equations" 

3 x/0j^ = 3f/30j^ + F(t) 3 x/30j<. (,20) 


where 


F(t) 


3F(t)/3 X 


A 

X 


The initial conditions of 3x/0j^ are zero except when 

A 

x(0) is identified as part of the parameter vector 0. 
In this case the initial partials have the value one, 

5 The li)<;elihood function for zero system noise is 

N 

J = log pCy . . y_/0) = - (1/2) S v: v 

j. n “t ^-2 

j=i j ^ _ 

Determine now the gradient of this function with 
respect to 0 

N T 1 

3J/30 = Z v: R~ 3V./30 (22) 

3=1 ^ ^ 
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where 3v./S© = - H 9 x./30 

6 Compute the information or sensitivity matrix 

N 

M :: 3^J/3Q^ = S Ov./ao)"^ 3V./30 

j=l ^ ^ 

The inverse M ^ of the information matrix 
provides a lower bound for the covariance of 
the updated parameter estimates. 

7 The updated parameter estimate is 

A 

0 = G + A 0 
o 

where A0 = (3J/30)^ 

8 Go now back to equation 18 with the updated 
parameter estimate and repeat the steps to 
equation 25, Reiterate until convergence of 
the information matrix and of the parameter 
vector is obtained, 

SINGLE BLADE PARAMETER IDENTIFICATION 

We first present a case where the sequential 
estimate is compared to the global estimate of the 
parameter y (blade Lock number). The excitation of 
the progressing mode is given by 


(23) 


(24) 


(25) 
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0 = 1,5 cos (o) H* l)t 

0 for t 4 12 

-( ,1/Tr)(t-12) for t > 12 

This excitation is somewhat different from that defined 
by equation 2 and illustrated in Figure 1, While 
equation 2 defines a constantly accelerated progressive 
excitation, equation 26 deviates somewhat from a 
constant acceleration. The difference is, however, not 
essential. The collective pitch angle is 0^ - 0, The 
rotor advance ratio is y = ,4, The blade frequency in 
the rotating system is P = 1.2. System and measurement 
equations for the single blade case are given by 
equations 14 and 15 of reference 1, 

The sequential estimate of y and of its covariance, 
was found from equations 10 and 11, using as initial 
estimate at time t = 12 the value zero, and using the 
initial covariance = 2000. The global estimate 
was determined from equation 12 setting the initial 
covariance to infinite, or = 0. Figure 3 gives the 

A 

estimate y and its covariance from the sequential method 
for the time t = 12 to 24 and also the value of the 

A 

global estimate, y = 5,05. The true value of this 
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parameter* is 5,00, From the covariance plot it is seen 
that at t = 16 the covariance is quite small as 
compared to the Initial value so that any further data 
processing will not appreciably change the estimate. 

A 

This is born out by the y plot. The global estimate 
for the time t = 24 agrees with the sequential 
estimate at this time, as it should be, since there 
is little difference between assuming the initial 
covariance as = 2000 or Pq = The method used 

for figure 3 is the same as that used in reference 1, 
The flapping response data were noise polluted by 
white noise with a standard deviation of = ,2, 

This is a large noise since the maximum flapping 
excursion is only 1.2. The polluted flapping 
deflections were then passed through a digital filter 
with the cut-off frequencies of 2,5 to 2.9. Rates 
of deflections and accelerations were determined from 
a Kalman filter with the erroneous value of y = ^.0. 
Sequential and global estimates were determined with 
the flapping deflections from the digital filter output 
and with the flapping rates and accelerations from 
the Kalman filter output. It is seen that this method 
worked very well in this case. The effect of the 
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digital filter cut-off frequencies is, however, 
substantial. The global method gave for the same 
case the following estimates 


Cut-off 

Frequency 

A 

Estimate y for t = 24 

b> 

C 

(0 

t 


1,7 

1.9 

5.27 

2.5 

2.9 

5.05 

3.5 

3.9 

4.89 


Because of the sensitivity of the estimate with 
respect to the cut-off frequencies of the digital 
filter it was decided to omit the digital filter and 
to use only the Kalman filter, however in an iterative 
way. 

The following example is again for a rotor advance 
ratio of y = ,4 blade frequency in the rotating system 
( 0 ^ = 1,2 and for a collective pitch setting of 9^ = 2°, 
Con.stant acceleration of progressive pitch stirring 
according to equation 2 is assumed, whereby t^ = 12, 

The polluted flapping response was processed with a 
Kalman filter using the initial parameter values of 
Y = 4,0 and 6 = yO^ = 8,0 as compared to their true 
values of 5.0 and 10,0 respectively. The two parameters 
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assumed to be unknown, y and 6, were determined by 
the global estimate equation 12 with P”^(0) = 0 and 
R = I. For the initial conditions of the Kalman 
filter B(0) was taken from the simulated measurement 
and 0(0) was taken with a 20% error. The following 
table gives in addition to the parameter values for 
4 iterations also the diagonal terms of the 
covariance matrix. There is some overshoot in both 
parameters and the convergence is not very good. 

Note that the covariance does not properly reflect 
the actual errors in the parameters.’ 


= .2, R = I 


Parameter 

^ 1 

5 = Y0o 


"6 

True value 

5.00 ! 

10.00 



Initial estimate 

4.00 

8.00 



Iteration 1 

(1.6 CPU 

5,66 

8.95 

.092 

.113 

seconds ner 2 
iteration ) 

4.76 

10.31 

.082 

.111 

3 

4, 52 

10 . 09 

.073 

.111 

4 

4,59 

9.93 

.071 

.111 


We now treat the same problem with the maximum 
likelihood method, and first assume the correct initial 
conditions 0(12) = -.921 and 0(12) = .986. The first 3 
iterations are given by 
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- .20 


Parameter 

Y 

5 

a J/Oy 

3J/a6 

True value 

5 . 00 

10.00 



Initial estimate 

4.00 

8.00 



Iteration 1 

(2.6 CPU 

4 . 9it 

9.70 

-6.11 

-5.84 

seconds per 2 

iteration ) 

4.91 

9 . 82 

.10 

-.34 

3 

4. 91 

9.82 

-.002 

, 001 


Though the computer' CPU tine is now per iteration higher 
than for the preceding method, the convergence is very 
rapid and the accuracy is very good. The Cram^-Rao 
lower bound for the parameter covariance matrix for 
the second and third iteration is 


M 


-1 


, 211 
-.036 


-.036 
. 365 


R = I 


In order to obtain the covariances for the actual 
simulated measurement noise = .20 one must multiply 
the above values by ,04, Thus the lower bounds for 
the parameter standard deviations are = .092, 

Cg = .12 which is quite close to the errors in the 


second iteration 
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In a second version of the same problem it is 
assumed that the initial conditions are unknown and 
must be included in the parameter identification. 
The following table gives the results of the first 
3 iterations. 




Parameter 

Y 

6 

8(12) 

8(12) 

True value 

5.00 

1 

10.00 

-.921 

.986 

Initial estimate 

4.0 0 

8.00 

-.750 

1.200 

Iteration 1 

4.96 

9.68 

-.89 

1. 02 

(4.3 CPU 





seconds per 2 

4.90 

9.85 

-.89 

1.01 

iteration) / 




** 

3 

4.90 

9.85 

-.89 

1.01 


The par>a*net d” covar’ance lower bounds are 




The computer CPU tine per iteration is once more increased, 
but the convergence is excellent and the accuracy of the 
estimate is very good. Again the inverted information 
matrix gives physically meaningful covariances for the 
unknown parameters. 


One can then conclude that for single blade 
identification with two unknown parameters y and 
5 = yQq iterative equation error estimation with 

updated Kalman filter gives the lowest CPU time per 
iteration, however because of the much faster convergence 
of the maximum likelihood method and because of its 
greater versatility (inclusion of initial response 
values in the set of unknown parameters), it definitely 
is preferable. For the following cases the maximum 
likelihood method will be used. 

MULTIBLADE PARAMETER IDENTIFICATION 

The preceding example is now treated with a multi- 
blade representation. The excitation in multiblade 
coordinates is given by equations 3 to 5 with to = 12. 

The pitch input is shown in Figure 2 assuming a 
constant progressing stirring acceleration given by 
equation 6, The multiblade responses ^o> ^ii are 

polluted from t = 12 to t = 24 with aero mean Gaussian 
noise of standard deviation .2, 

We first identify the parameters Y and 6 s yQq 

• • • 

with given initial values of 3 q » 3o » 3is 

at t = 12, The responses are obtained with the periodic 
system equations. Hov:ever, for the parameter identi- 
fication, constant coefficients are assumed so that 
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modeling errors are present. The initial states at 


t = 12 

are: 





^ = 

.880, 

Bj 

= -.983, Bjj 

= 1.698, 

B, = .060 

d 

• 


• 

• 


A 

6 = 
o 

.038, 


= .034, Bjj 

= .03 , 

B, = -.030 
d 

refers to 

the 

reactionless 

mode of a 

four-bladed 

rotor , 

Bo is 

the 

coning angle 

, Bj represents forward 


tilt and left tilt. For the constant coefficient 

approximation 3^ does not couple with the other 
states and can be omitted. The multiblade equations 
for a four-bladed rotor are given in reference 9, 

The exact and the noise polluted responses Bj, 3^^ 
are shown in figure 4, 

Using the above initial conditions one obtains with 
the maximum likelihood method the following values 


Parameter 

Y 

6 

True Value 


5.00 

10,00 

Initial Estimate 

i 

2.50 

5.00 

Iteration 

1 

4.97 

8.46 

(4 CPU seconds 

2 

5.00 

9.62 

per iteration) 

3 

i 

4.99 

9.65 
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The lower bound for the parameter covariance in the 
3rd iteration is 


M 


-1 


.125 

.003 


" 

.003 

.269 


, R = I 


The result as compared to the preceding single blade 
identification is better for y and worse for 5. Note 
the assumed large errors of 50% in the initial para- 
meter estimate. 

In order to assess the effect of errors in the 
initial states at t = 12, it is now assumed that all 
initial states are zero. One then obtains: 


3o = Bq = = Bjj = Bjj = 0 at t = 12 


Parameter 

y 

6 

True value 

5.00 

, 10.00 

Initial estimate 

2.50 

5.00 

It era tion 1 

5.68 

3.U-5 

(4 CPU seconds 
per iteration) 2 

4 . 87 

10.27 

3 

4.92 

10.26 
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For R = I the inverted information matrix is 

ol51 ,043 

.043 .278 

The result shows that large errors in the initial states 
can be tolerated. 

Finally the 3 initial deflections have been assumed 
as unknown parameters, increasing the number of para- 
meters to be identified from 2 to 5, 'The initial 
values for the response rates were assumed to be zero: 

8 = 8-r = = 0 at t = 12. The initial values of 

8^, 8j, are given in the table: 


Parameter 

Y 

1 5 

1 

8^(12: 

3i(12^ 

6 ji(12 

True values 

5.00 

10.00 

.880 

-.983 

1.698 

Initial estimates 

2.50 

^ r\ r* 1 

1.00 

-1. 00 

2.07 

Iteration 1 

(6 CPU seconds 

4.94 

8.25 

.979 

-1.039 

1. S26 

per iteration) 2 

5.01 

9.97 

.920 

-.985 

1.732 

3 

5.00 

9.93 

,920 

-.984 

1. 730 
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The inverted information matrix for y and 5 is: 



It is seen that now the convergence is much faster, 
the second iteration being even better than previously 
the 3rd iteration. However, the CPU time per 
iteration is increased from 4 seconds to 6 seconds. 

The following table compares the results of the 
various methods. The last 4 rows refer to the maximum 
likelihood method. 



The number of iterations indicated in the table is that 
for which convergence has been achieved. The iterated 
equation error estimation with updated Kalman filter 
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needs the lowest total computer effort, however, the 
accuracy of the estimate is worst for y. The maximuni 
likelihood estimation, due to faster convergence, needs 
only moderately more computer effort and yields better 
accuracy. For the single blade the identification of 
the two initial states, 3(12) and 6 ( 12 ), results in the 
same good accuracy of the parameters as when using 
the correct initial conditions. 

The multiblado identification - despite modeling 
errors by omitting all periodic terms in the equations 
of motion and despite assuming zero initial con- 
ditions - also converges rapidly and provides only 
slightly less accuracy of the parameters. The multi- 
blade analysis with identification of the 3 initial 
defle-ctions gives the best accuracies, however with 
more than twice the computer effort. Note that the 
convariance estimate of the Kalman filter method is an 
order of magnitude greater than for the maximum 
likelihood method. The latter values have been obtained 
from M ^ for R = I by multiplication with ,04, since 


R 


, 04 
0 


0 

,04 


They represent lower bounds of the 


actual parameter covariances but agree fairly well with 
the errors found for the identified parameters. This 
is not true for the covariance estimates of the Kalman 
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filter method. The latter method was attempted 
also for the multiblade representation but was 
found to require impract ically long computer runs 
per iteration combined with a lack of convergence. 

The multiblade identification studies were, therefore, 
limited to the maximum likelihood method. 

DYNAMIC WAKE PARAMETER IDENTIFICATION FOR ZERO ADVANCE 
For the case of zero advance ratio reference 10 
gives a mathematical model of the dynamic rotor 
wake that has 2 constants: The quasi~steady wake 

number L that relates the pitching and rolling 
moments of the rotor to the sine and cosine components 
of the wake velocity, and a time constant t. The 
third rotor constant is 

A = Y/8 

where B is the tip loss factor and y the blade Lock 
number. The' inflow angle from the rotor wake is 
assumed constant over the radius and represented 
for the kth blade by 

X, = X + X_ cos + X^_ sin jj;, 
k o I ^k II k 

In a linearized analysis for zero advance ratio the 
coning mode is uncoupled from the tilting modes. The 
system equations of reference 10 are in matrix notation 


RATIO 


(27) 


( 28 ) 
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There are 6 state variables: 8^ , 8^ , 8jj , X^j. , 

The flapping angles of the 4 blades 8^^^, 82* 3g, 8^^ hre 
measured in the rotating system and combined into 
^l”^3 measurement equations - that is 

the relation between the measured quantities and the 
system state variables - are given by 



where p , rj are measurement noise of 8 , and 8 « 
12 ml T m2 

respectively. Thus only 2 of the 6 state variables 


( 30 ) 
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are measured. The parameters to be identified are 
A, AL/t, 1/t, and the blade natural frequency in 
the rotating system is assumed to be given: 
to 2 ^ — 1,2, 

The pitch stirring input is again described by 

equations 4 to B (progressing mode). The rotor 

wake causes an increased stabilization time. 

Therefore t =70 was selected instead of t = 12, 
o o 

The measurement noise assumed to be vrhite 

and Gaussian with the standard deviations 


a 


ml 



= ,10 


The reduction to one half the noise amplitude as 

compared to the preceding examples is justified by 

the fact that 6 , and 6 . are substantially smaller 
ml m2 

than figures 4 and 5,* The true values 

of the 3 parameters for which the response was 
determined are 

A = .500 AL/t = .250 1/x = .125 

These values are in the range expected to be found 
for the model rotor. They will depend essentially 
on the collective pitch setting. In preliminary 


3m 


*Note that figure 5 refers to a 


.05 
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computer runs it was found that the initial estimate 
of 1 /t should not be zero. 

Using the true initial conditions one obtains 
with the maximum likelihood method 


Paranet cr 

A 

1 

AL/t j 

1 1/t 

True value 

.500 

.250 

! .125 

! 

Initial estimate 

. aoo 

. 200 

1 

. 250 

1 

Iteration 1 

.500 

.290 I 

. 086 

(3.8 CPU 


1 

1 

seconds oer 2 

. 48 0 

. 253 ! 

1 .125 

iteration ) 



1 

3 

.487 

. 255 ! 

j 

! . 128 

1 


4 .487 .255 ‘ .128 


With the inverted information matrix for the 4th 
iteration 


052 

.059 

n 1 ^ ^ 

• J ^ o 

j 

0 59 

.094 

. 023 1 
1 

013 

.023 

. 010 J 


R = I 


Despite the fact that only 2 out of 6 state variables 
are measured, the parameter identification is very 
good and the second iteration has almost converged. 

The following case is the same as before except 
for regressing excitation (+o)) instead of progressing 
excitation (-w) 
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Parameter 


! 

! 

A 

AL/t 

1/T 

True Value 



o500 

.250 

.125 

Initial Estimate 



.400 

.200 

.250 

Iteration 

1 


.524 

.262 

.087 

(3.8 CPU 

2 


. 540 : 

.264 

.134 

seconds per 

3 


. 536 j 

.264 

; ,139 

iteration ) 

4 


. 536 1 

,264 

' .139 

with 







• 

167 

.179 

.029 ' 


= 

• 

179 

.256 

. 045 

II 

M 


• 

029 

. 045 

.013 _ 



Though convergence is good and the second iteration 
has almost converged, the errors are now larger 
than before as expressed also by the M ^ matrix. 

The physical reason why regressing modes are less 
suited for rotor wake identification is the fact 
that at (0 = ,2 the excitation is in resonance with 
the regressing flapping mode. At this condition 
no dynamic rotor wake exists since aerodynamic 
excitation and aerodynamic damping cancel each other. 
Since regressing mode transients include a frequency 
region with a weak dynamic rotor wake, the identi- 
fication of the wake parameters is less good than for 
progressing mode transients. 


airpSODUCIBILinr OF THS 
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Finally the same case is treated with very large 
errors in the initial parameter estimates. Now the 
4 initial displacements are identified in addition 
to the 3 parameters, resulting in 7 instead of 3 
unknown parameters. The initial rates are assumed 
to be zero: = X^^. = 0 at t = 70, 

Using progressing excitation one obtains 


Parameter 

A 

AL/t 

1/t 

Bi(70) 

«„(70) 

Xj(70) 

Xjj(70) 

True Value 

. 500 

.250 

.125 

.497 

.188 

-.874 

-.331 

Initial Estimate 

.400 

. 067 

. 083 

. 601 

.305 

-.531 

-.268 - 

Iteration 1 

,443 

.184 

,150 

.500 

.211 

-.864 

-.309 

(6,3 CPU seconds 2 

.471 

.233 

.115 

.497 

.210 

-.859 

-.322 

per iteration) ^ 

.474 

,235 

.113 

.497 

.211 

-.870 

-.327 

4 

. -1 , . 

! .474 

1 

.236 

,120 

.497 

.211 

-.871 

-.327 

!I diagonal 4th 

iteration, R = I 

^ .121 

.184 

.055 

. 088 

. 095 

1.069 

.796 


The convergence is again good and the second iteration 
is almost converged. The errors in the 3 parameters 
are somewhat larger than before but still acceptable. 
The errors in the initial conditions, except for 3jj, 
are also small, though the H ^ matrix shows larger 



covariances for and The initial estimate of 

the initial displacements 3 3j j, X j, X at t = 70 
are not arbitrary but are determined from the 
equilibrium equations before the beginning of pitch 
sTirring, using the initial estimates of A» AL/t 
and I/t. Since the true initial conditions for Xj 
and Xjj are unknown, their inclusion in the identi- 
fication is necessary. and are measured, thus 

their initial conditions are known within the accuracy 
of the measurement and their inclusion in the 
identification may not be required if the measuring 
error is small. The computer time for identifying 
all 4 initial displacements (3p X^., X^^ at t = 70) 

is hoxirever quite moderate, so that their identi- 
fication is practical. 

On the basis of the computer experiments the 
method to be used for the rotor model tests at 
zero advance ratio is tlien as follows: 

1. Estimate the parameters A, AL/t, 1 /t 

2. Determine from the steady state equations 29 

with zero rates and zero accelerations the 
initial values for time t = t^, 

using the estimated parameters. 

3. Assume zero initial rates: 




0 at t 
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Apply the maxitnuni likelihood algorithm to 
equation 29 and determine iteratively the 7 
parameters A, AL/t, 1/t, l^(tQ), 

Xll(to). 

SIMULTANEOUS MEASUREMENT AND SYSTEM NOISE 


Though the simplified maximum likelihood 
algorithm as currently used in aircraft parameter 
identification work and as defined here by equations 
18 to 25 does not provide for system noise, such 
noise is unavoidable in full scale and wind tunnel 
tests. The origin of the system noise is either 
in atmospheric or wind tunnel flow turbulence, or 
in modeling errors. In order to assess the detri~ 
mental effects of the system noise on the quality 
of the maximum likelihood estimates, the preceding 
example was recomputed under the assumption that 
equation 29 has an additional noise term on the 
right hand side of 


A 

0 


0 

A 



-(AL + 1)/t 0 

0 -(AL + 1)/t 



(29a) 


This means that we have noisy inflow components 
For a zero mean computer generated Gaussian 
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sequence with standard deviations ” ^^2 ~ 

assumed. Figure 6 shows the noise polluted inflow 
components 1^, addition noise was added 

as before in the measurement equations 30, however 
now with standard deviations O- - = Oo„_ = .05 which 
is one half the previously assumed value and which is 
more representative of the expected measurement errors. 

We first determine the effect of the lower 
measurement noise on the parameter identification: 


®3ml ' '^3m2 




Parameter 

A 

AL/t 

1/t 

6jC70) 

6jj(70) 

lj(70) 

Ajj(70) 

True Value 

. 500 

.250 

.125 

.497 

.188 

-.874 

-.331 

Initial Estimate 

.400 

.067 

. 083 

.601 

.305 

-.531 

-.268 

Iteration 1 

.455 

.189 

.159 

1 , 504 

.195 

-.880 

-.330 

2 

.483 

.238 

.119 

. 500 

.194 

-. 867 

1 

-.339 

3 

.485 

.239 

.123 

, 501 

.194 

-.874 

-.341 

4 

.485 

.239 

.123 

i 

.501 

.194 

1 

* -.874 

-.341 


It is seen that the accuracy of the parameters is 

improved as compared to the preceding case of 

drt 1 = cr„ _ = .10, as it should be. 

3ml 3m2 
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Combining now the measurement and system noise 
one obtains: 


'^3ml " '^3m2 



.10 


Parameter 

A 

AL/t 

1/x 

Bj(70) 

B„(70) 

Xj(70) 

Xjj(70) 

True Value 

. 500 

.250 

.125 

.497 

1 

.188 

-.874 

-.331 

Initial Estimate 

.400 

.067 

.083 

. 601 

. 305 

-.531 

-.268 

Iteration 1 

.455 

.189 

.153 

. 505 

.188 

-.860 

-.319 

2 

.482 

.236 

.115 

.502 

.188 

-.847 

-, 324 

3 

.483 

.235 

.118 

. 502 

.190 

-.853 

-.325 

4 

.483 

.236 

.118 

. 502 

.190 

-.854 

-.325 


It is seen that the system noise indicated in figure 6 
has only a small detrimental effect on the accuracy 
of the identified parameters. Of all computer 
experiments of parameter identification we have shown 
here only the results of one computer run. Since these 
results are dependent on the computer generated random 
sequences, we are actually dealing with identified 
parameters that are random variables. A repeat run 
will, however show only minor differences in the 
identified parameters, so that the presented results 
in comparison to the true values can be considered 
as typical. The parameter that suffers most in its 
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accuracy when system noise is added is the reciprocal 
time constant 1/x. Even so the error is only 5,5% 
which would appear quite acceptable for all 
practical purposes, 

CONCLUSIONS 

The computer experiments presented here can be 
used to establish preferences with respect to various 
alternatives in the rotor parameter identification 
schemes . 

1 . Iterated Equation Error Estimation with Updated 
Kalman Filter vs. Maximum Likelihood Method 
For single blade parameter identification from 
pitch stirring transients the equation error 
method applied in an iterative form using a 
Kalman filter with the latest parameter updates 
worked well and required the least computer CPU 
time. For raultiblade parameter identification 
this method became impractical because of slow 
convergence and high computer CPU time. The 
maximum likelihood method worked well both for 
single blade and multiblade applications, though in 
case of single blade identification it requires 
somewhat more computer CPU time. The parameter 
covariances from the maximum likelihood method are 
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clearly superior to and more meaningful than the 
covariances determined with the equation error 
method. The maximum likelihood method also gave 
good parameter identifications in the presence 
of both measurement and system noise, though 
most of the computer experiments were conducted 
with measurement noise only. Overall one can 
conclude that the maximum likelihood method in 
its simplified form in which system noise is 
not modeled, is for the applications studied 
superior to the equation error method and thus 
will represent the method of choice for the 
parameter identification from wind tunnel rotor 
model tests. 

2. Fixed vs. Identified Initial State 

Since prior to the initiation o'f pitch stirring 
the rotor is approximately in a steady condition 
as far as the multiblade coordinates are concerned, 
one can always assume their initial rates as zero. 
For two parameter identification in forward flight 
(u = ,4) the use of zero initial values for all 
states gave good identification results. Including 
the initial displacements in the identification - 
while retaining initial zero rates of displacement - 
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greatly improved the convergence and the 
accuracy, though at a cost of about twice the 
computer CPU time per iteration. For the mul- 
tiblade analysis it appears nevertheless 

practical to identify the initial displacements 

\ 

in order to improve convergence and accuracy 
of the parameters. For the single blade analysis 
the initial value of the flapping rate is not 
zero and must be identified, 

3 , Single Blade vs. Multiblade Analysis 

If rotor wake lag is included, only a multiblade 
analysis is possible. Without rotor wake lag and 
using an "equivalent" Lock number to describe 
the wake effects, a single blade analysis is 
possible. Initial flapping deflection and rate 
of deflection should be identified. The multi- 
blade analysis with zero initial rates and 
identified initial deflections is more accurate 
though more demanding in computer CPU time per 
iteration. The multiblade analysis is thus 
preferable also in those cases where a single 
blade analysis is possible, 

4 . Progressing vs. Regressing Pitch Stirring 
Progressing pitch stirring leads to more accurate 
rotor v/ake parameter identification than regressing 
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pitch stirring, because the regressing stirring 
transient passes through the flapping resonance 
at which no dynamic rotor wake exists. 
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Chapter 2 
FIGURE CAPTIONS 


Figure 1 
Figure 2 
Figure 3 


Figure 4 


Figure 5 


Figure 6 


Tine History of Blade Pitch, Rotating 
Reference System, Progressing Mode 

Time History of Blade Pitch, Multiblade 
Representation, Progressing Mode 

Estimate y and its Variance for 
P(0) = 2000, y(0) = = *2, 

0Q = 0, u = .4, coi = 1.2, Progressing 
Excitation 

Exact and Noise Polluted Multiblade 
Responses (Jjj- to Inputs Shown 

in Figure 2 for y = ,4, 0 = 2°, 

Y ~ 5.00, toj_ — 1.20, *" ~ ^Bii *" 

Exact and Noise Polluted Multiblade 

Responses S , , 6 ^ to Inputs Shown in 
^ ml m2 

Figure 2 for y = 0, A = .500, 

AL/t = , 250, 1/t = .125, = 1.20, 

'^Bml ^ *^6m2 ~ 

Exact and Noise Polluted Inflow Components 
Xl, Xjj (system noise) for y = 0, 

A = .500, AL/t = .250, 1/t = .125, 

till = 1.20, a, = a, = .10 

^I ^II 
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CHAPTER 3 


OPTIMAL DATA UTILIZATION FOR PARAMETER IDENTIFICATION 
PROBLEMS WITH APPLICATIOH TO LIFTING ROTORS 


K, H, Hohenemser and D, Banerjee 


ABSTRACT 

A method is developed for optimal data utilization 
in the maximum likelihood identification of systems 
without process noise, based on stipulated upper bounds 
of parameter covariances. The method is applied to a 
case of simulated transient wind tunnel testing of 
lifting rotors. 
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INTRODUCTION 

Methods for state and parameter estimation from 

12 3 4 

transients are widely used in aircraft testing. » » * 

The problem is to obtain optimum estimates (based on 
certain performance criteria) of initial states and of 
unknown parameters (derivatives) from noisy measurements 
of some inputs and response variables. In most cases 
of airplane parameter identification a constant system 
is used as an analytical model. For lifting rotor appli- 

5 

cations a periodic system model is required. State 
and parameter identification from transients looks 
promising also for wind tunnel testing and may well 
drastically reduce the amount of test efforts as compared 
to e.g,, frequency response testing of lifting rotors. 

The following study was performed in preparation for 
transient wind tunnel tests with a lifting rotor model. 

In the most general case the identification method 
has to account for the following types of uncertainties: 
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(a) Modeling errors which originate from the difference 
between the mathematical model and the actual 
physical phenomenon, 

(b) Uncertainties from input noise, e,g,, turbulence 
(in the atmosphere or in the wind tunnel), 

(c) Uncertainties from measurement inaccuracies. This 
includes measurement bias and noise and also incom- 
plete measurements when some of the input and 
output data are missing. 

Three basic methods have been applied to the air- 
craft identification problem. If the measurements are 
a linear function of the unknown parameters, the 
classical least squares regression or equation error 
method is applicable. This method was widely used in 
the early stages of aircraft parameter identification. 
For the application of this method the time histories of 
all the states and their derivatives are needed together 
with the input variables. Since the equation error 
method becomes less reliable the noisier the measure- 
ments are, it has been replaced in the last decade by 
less restrictive methods that allow considerable 
measurement noise and that work also when some states 
have not been measured. One of the most widely used is 
the maximum likelihood method, which, in the absence of 
modeling errors and input noise, reduces to what is also 
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12 4 

called the output error method. * * Another way of 
overcoming the limitations of the least squares 
•regression method is the application of the extended 

3 

Kalman filter. This note is applicable to data 
analysis with the output error method. This method is 
not necessarily limited to smooth air testing. Atmos- 
pheric or wind tunnel turbulence may be included if 

6 7 

appropriate input measurements are taken. * 

In aircraft or wind tunnel transient testing the 
question comes up as to what kind of transient should 
be selected. If the transient is too short, the 
parameters will be identified with inadequate accuracy. 
If the transient is too long, an unnecessary amount of 
data must be processed. The Question we pose here for 
the maximum likelihood method is: Given a reauired 

accuracy of the parameter estimate, and given an input 
function, what is the minimal quantity of measured data 
necessary to achieve this accuracy? From a survey of 
the pertinent literature it appears that this question 
has not been asked before, much less a solution 
presented. However, there are some recent studies 
where certain criteria were used to define an optimum 
input. We will first briefly discuss two of these 
optimal input proposals, and then proceed to develop 
the method of optimal data utilization for a given type 
of input. 
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TWO PROPOSALS FOR OPTIMAL INPUT DESIGN 

General questions of input design are: 

(a) What type of input function should be used? 

(b) For what time period should the response data be 
processed to enable identification of the system 
parameters with a specified accuracy? Are certain 
time periods of the response particularly rich in 
information contents and should they, therefore, be 
preferably used? 

There usually are some constraints on the input design 
like amplitude constraints, smoothness constraints 
(step or impulse inputs are mathematical idealizations 
but often practically not realizable), instrumentation 
constraints, and constraints imposed by the selected 
analytical model that usually filters out the higher 
frequency contents of the input. 

Analytical solutions of the problem of optimal 
input design require the minimization of a cost function. 
Stepner and Mehra^ use the sensitivity of the system 
response to the unknown parameters as the performance 
criterion for optimal input design. The time of the 
transient is assumed to be fixed. Thus questions (b) 
are not involved. The measurement equation is 

y^(t) = y(x, 0, u, t) + v(t) (1) 

X is the state vector, 0 the parameter vector, u the 
input vector, v(t) the additive measurement noise with 



covariance matrix R, We write the Taylor expansion 
- with respect to the parameter 0 about the a priori 
estimate 0^ of 0 and neglect higher order terms: 


So* y(x, 0 ^, u, t)(0-0o) 

+ v(t ) 

In the output error method (0 - 0 q) is determined by a 
least squares solution of Eq. (2) for a fixed time 
period (tg* tf)» for a high degree of accuracy in 
determining (0 * 0g) the sensitivity function 3y/30 
must be large. The scalar performance index selected 
in Reference 1 is 


( 2 ) 


J = Trace (WH) 


( 3 ) 


where 


M = 


f 

I ' 


ay/ae)"^ R“^(ay/8a)dt 


( 4 ) 


Due to the introduction of in M the performance 

criterion favors the measurements which are more accurate. 
The weighting matrix W is based on the relative impor- 
tance of the parameter accuracies. 

Assuming now linear system and measurement equations 


x(t) = F x(t) + G u(t) 

y (t) = H x(t) + v(t ) 
m 


(5) 

( 6 ) 
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together with an “energy constraint” for the input 

E = / u dt (7) 

the optimum input u can be determined as a two point 
boundary value problem whereby the Hamiltonian includes 
the term 

(l/2)Po (u*^ u - E/tf) . 

Vo is the time invariant Langrange factor (scalar) to 
be evaluated from the Euler differential equations of 
the optimization problem. It should be noted that the 
"energy constraint" Eq, (7) has no physical signifi- 
cance but is a convenient device to obtain smooth 
input functions. Physically the input will usually 
be limited by amplitude rather than by the quadratic 
criterion (7) and quite different "optimal" inputs can 
then be expected. 

Chen® attacks the problem of optimal input design 
in an entirely different way as a time-optimal control 
problem by minimizing 


"^-f 



under the following constraints; 





System equations 

X = f(x s u, e s t) , x(t ) = X (9) 

o o 

Sensitivity equations 

3A/3Xq = Of/3x)Ox/3xo) , Dx(tQ)/3Xo = I (10) 

3x/30 = (3f/3x)(3x/36) + 3f/30, 3x(tQ)/30 = 0 (11) 

t 

Information matrix equations 

= ~H"^(3v/3e)^ r“^(3v/30 )h”'^ (12) 

where V is the innovation; 

V = yjn - y (13) 

and where the information matrix H is given by Eq. 4, 

Finally Chen assumes an amplitude constraint 

luj < V (14) 


and he prescribes the trace of the information matrix 
for time tf 

C..(t^r) = a? (15) 

One can show that for linear input u(t) into the system 
equation and for an input matrix independent of any 
unknown parameter, the optimal input is of the *'bang~bang” 
form between the amplitude constraints. The solution of 
this problem requires a computer search which was not 
performed in Reference 8, Rather an arbitrary set of 
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bang-bang inputs in the form of Walsh functions was 
shown to result in a specific case in lower values of 
M“^(t^) (given tf) than those obtained by using 
Mehra’s ’’optimal input”. This apparent contradiction 
can be explained by the differential equation (12) 
governing M~ . For a particular value of the 

rate of decrease of with time is dependent on 

all elements of 

Oy/30)’’'' R”^Oy/30) 

while Hehra, in his criterion (3), optimizes only the 
trace of WM, 

While the input amplitude constraint Eq, (14) 
used by Chen is physically more significant than the 
quadratic constraint Eq, (7) used by Hehra, the actual 
constraints are usually still more complex. In cases 
of airplanes or lifting rotors one usually wishes to 
limit the response to the linear sub stall regime, 
since the analytical model to be identified is often a 
linear one. The stall boundary is, however, a complex 
function of the input and cannot be represented by an 
amplitude constraint for the input transient. This is 
particularly relevant for the lifting rotor, so that 
neither the Mehra nor the Chen input optimization 
criteria are useful for lifting rotor applications. 
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quite apart from the excessive computer effort involved 
in obtaining the optimal inputs. Furthermore, the 
input matrix usually contains unknown parameters. In 
this case Chen's optimum solution would not be of 
the bang-bang type and would be still more difficult 
to obtain. For all of these reasons it was concluded 
that at the present state of optimal input design 
methods an attempt to compare our selected inputs with 

an "optimum input" would not be practical. Instead, a 
more limited approach has been taken described in the 

following section, 

OPTIMAL DATA UTILIZATION FOR GIVEN INPUT FUNCTION 
We first point out the difference between the 
continuous and the discrete case. In the maximum 
likelihood method (output error method for zero process 
noise) using the Newton-Raphson approach with quasi 
linearization, one obtains for the parameter update 
increment the following expressions: 

Continuous case; 


A6 = 


/lv,T 




L-t, 


f (IQ) 


Discrete case; 


A0 “= 


r N 
S 

i=l 




-1 / 5 V 


1 -1 


^30 1 


N 

E 

i=l 


,3v.T 


EEPEODXICIBILiry OF M 
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(H) 

^30 ■' 


( 17 ) 
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The Cramer-Rao lower bound has been defined only for 

a vector of sampled measurements and not for the 

2 9 

continuous case, * 

For high sampling rate one can define an approxi- 
mate differential equation for M from Eq, (17) in the 
following way: 


Set S. A (3u/36 ) . 

1 "• 


(18) 


then 


M 


N 

= E 
i=l 




N 

S 

i=l 



(19) 


As N increases. At gets smaller and the right hand side 
of Eq. (19) can be approximated by 


M s (1/At) 



= NAt 


T -1 
S R 


S dt 


(20 ) 


s At 


■*^f 

I 

*'t 

t- '-o 


T -1 

S R S dt 


-1 


-1 


Taking the derivative of M with respect to t^ 


dM"^/dt^ = -M'^OM/atf )m"^ 


( 21 ) 


or with Eq. (20): 

dM"^/dt^ = -d/At)M~^ S M 


( 22 ) 
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The point is that even in a continuous formulation the 
tine increment At between samplings must occur. 

Eq, (21) is the correct formulation for the Cramer-Rao 
lower bound of the covariance matrix for the parameters. 

We can now use the approximately valid differential 
equation (22) to obtain some insight into ways of best 
data utilization. Let us assume that we wish to pre- 
scribe certain values for the parameter standard 
deviations cj£ and that we wish to compare the Cramer-Rao 
lower bound with these standard deviations. Since we 
are dealing not with the unknown actual parameter 
covariances but only with their lower bounds, we should 
apply some conservatism to the selected aj, that is we 
should select smaller than we really need for the 
specific data processing case. We thus require 

0 < Mtf (i, i) £ a? (23) 

whereby is the value of at time t^. For 

non-zero values of S, the right hand side of Eq, (22) is 
negative definite and hence (i, i) are monotonically 

decreasing functions of tf. There will thus be a 
minimum time for which the constraints of Eq. (23) are 
satisfied. 

Another way of reducing the amount of measured 
data for the parameter identification is to select for 
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the data processing those time periods for which the 
components of the matrix 



have significant values. From Eq. (21) it follows 
that then the Cram^r-Rao lower bound will be 

particularly small. The components of M ^ also 
decrease with decreasing time element At between samples. 

Since it is impractical to use for the integration 
of Eq. (22) infinity as initial condition, it is 
recommended to determine M“i for a small time period, 
say for N = 10, from Eq. (21) and integrate Eq. (22) 
with the solution to Eq, (21) as initial conditions. 

Since S includes parameter estimates, one needs a 
preliminary estimation of the unknown parameters in 
order to use Eq, (22), 

APPLICATION TO A CASE OF LIFTING ROTOR PARAMETER IDENTIFICATION 
Using the simplest analytical model of a lifting 
rotor, a straight blade flapping about the rotor center, 
one has in a rotating frame of reference for the flapping 
angle g the following equation, 

3 + (Y/2)C(t)| + (y/ 2)[P^ t K(t)]g = (Y/2)[mQ 9 + A] 

(24) 

where: y is the blade Lock number 

P is the blade flapping natural frequency in 
the rotating system 
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6 is the instantaneous blade pitch angle 
X is the non-dimensional normal inflow 

One rotor revolution corresponds to t = , For 

neglected reversed flow effects, zero root cut-out 
and with tip loss factor B, the functions C(t), K(t), 
m 0 (t), in^(t) in terms of rotor advance ratio u are:^^ 

C(t) = (1/4 )b'^ + (1/3)B^ y sin(t) 

K(t) =(1/3)B^ y cos(t) + (l/4)B^ sin(2t) 
m^Ct) = (1/4)(B^ + s2y2) + sin(t) 

- (l/4)B^y^ cos ( 2 t) 

m;^(t) = (1/3)B^ + (1/2)B^ y sin(t) 

In the numerical analysis, we use B = 0,97, A 

simple Improvement of this analytical model that takes 

into account blade bending flexibility is possible. 

In transient conditions the inflow X includes the 

dynamic rotor wake in a complicated form. As a first 

approximation of dynamic rotor wake effects one can 

use, instead of the actual blade Lock number an equivalent 

12 

smaller value of Y» Such an approximation can be 

expected to be satisfactory if the transient is 

relatively slow. For transients with high frequency 

13 

contents, this approximation is invalid. 


(25) 

(26) 

(27) 

(28) 
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Due to rotor induced cross flow in a wind tunnel, 
the inflow parameter X will usually not be well known. 
In addition, the aerodynamic pitch angle 0^, due to 
airfoil inaccuracies and pitch setting errors, will 
also not be well known. For the wind tunnel tests 
considered here, we assume X = 0 and use the collective 
pitch setting 0^ as an unknown parameter to be 
determined from the blade flapping measurements. In 
addition we have a transient blade pitch input 6 
assumed to be known. The problem then is to determine 
from blade flapping transients caused by blade pitch 
inputs, the equivalent Lock number y and the 
equivalent collective pitch setting 0^, 

The lifting rotor wind tunnel model described in 
reference 14 allows excitation of progressing and 
regressing flapping modes at various frequencies. By 
a minor modification of this model, progressing or 
regressing transients can be excited. One can describe 
such inputs as pitch stirring transients. In a 
helicopter, this would amount to cyclic stick stirring, 
whereby the amplitude of the cyclic pitch would remain 
constant while the frequency of the stirring motion 
changes. The blade pitch input for such a stirring 
transient is selected to satisfy the equation: 


0 = 1.5 cos[&)(t - tg ) + t"] + 9o 


(29) 
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where oi is the angular stirring speed in the sense of blade 
rotation 

t =12 and 
o 

0 0 ^ t < to 

0) = 

-(0.1/7r)(t “ t^) to 1 t 1 T 

The meaning of these input equations is the fol- 
lowing- At time t = 0, a step lateral cyclic pitch 
input of 1,5 degrees is imposed. At time t = 12, 
the response to this input is approximately stabilized. 

At this time the pitch stirring acceleration of 
-(0,l/ir) is introduced which leads to 'a progressing 
flapping excitation. The identification starts at 
t = 12 with the pitch stirring transient. 

Test results with transient pitch stirring inputs 
will be presented at a later time when they become 
available. Here we are concerned with the problem of 
designing the tests in such a way that the test data 
will be sufficient to determine the two unknown para- 
meters Y and 0 q with good accuracy, i,e,, to determine 
a suitable value of T that allows an accurate identi- 
fication of parameters. 

The simulated identification analysis was performed 
under the assumption of a random zero mean white noise 
sequence superimposed on the analytical flapping transient. 
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This transient was obtained for 8^ = 2°, y = 0,4 and 

Y = 5.0, For convenience, the parameters 6 = Y^o 

Y instead of 9^ and y were identified. 

System and measurement equations corresponding to 
equations (5) and (6) are: 


• 

"^r 


0 

1 


-Xi 0 






+ 



_ -(y/2(P^ + K(t ) ) 

-(Y/2)C(t )_ 


-^2 (Y/2)me^_ 


ym 


[1 



+ v(t ) 


(30) 

(31) 


where 

E (v(t)} = 0 E[v(t) v(t)] = 0.2 5(t - t) 


and = [3 3] 

ANALYSIS OF RESULTS 

We first show in Table 1 the effect of data length 
on the narameters and their associated M“i(i,i) values. 
The iteration of the maximum likelihood method was 
begun with a 20 percent error in the parameter values. 
It is seen that a data length of t = 12 - 14 is quite 
inadequate, a data length of t = 12 - 18 gives 
reasonably good parameters, while a data length of 
t = 12 - 24 is much better and leads to very small 
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lower bounds of the parameter covariance matrix. 

Fig, 1 shows the correct flapping response together 
with the simulated measurement data. Pitch stirring 
is initiated at t = 12. Figs. 2 and 3 show 
and from Eq. (22) between t = 16 and t = 24. 

Two curves are plotted, one for the initial crude 
estimate of the parameters (y = 5 = 8 and one for 

the final estimate of the parameters for t = 24, 

(y = 4,91, 6 = 9,83), The two curves are in this 
case not much different. Mote the steep descent of 
the curves to about t = 17,5. It would, therefore, 
be not acceptable to use the data up to less than the 
time t = 17.5. However, there is another descent to 
t = 23,0, causing the improvement shown in Table 1. 
From Figs, 2 and 3 it is clear that the selection of 
T = 24,0 is a good one, that the use of fewer data 
would result in substantial decrease in parameter 
accuracy, and that the use of additional data is 


unnecessary 
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FIGURE CAPTIONS 

Fig, 1 Simulated measurements of 6, - ,2 

Fig. 2 Cramer-Rao Lower Bound M~^(y) vs. time 

Fig. 3 Cramer-Rao Lower Bound H"^(5) vs. time 

Table 1 Iteration vs. Time Variation of Parameters and 
Their Cramer-Rao Lower Bounds. 
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